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ABSTRACT 

Ground-based observatories have been collecting 0.2 — 20TeV gamma rays from blazars for about 
twenty years. These gamma rays can experience absorption along the line of sight due to interactions 
with the extragalactic background light (EBL). In this paper, we show that the gamma-ray optical 
depth can be reduced to the convolution product of an EBL kernel with the EBL intensity, assuming 
a particular form for the EBL evolution. We extract the absorption signal from the most extensive 
set of TeV spectra from blazars collected so far and unveil a broad-band EBL spectrum from mid¬ 
ultraviolet to far infrared. This spectrum is in good agreement with the accumulated emission of 
galaxies, constraining unresolved populations of sources. We propose a data-driven estimate of the 
Hubble constant based on the comparison of local and gamma-ray measurements of the EBL. After 
setting stringent upper-limits on the redshift of four TeV blazars, we investigate the 106 gamma- 
ray spectra in our sample and find no significant evidence for anomalies. The intrinsic TeV spectra 
are not harder than their GeV counterpart, and no spectral upturn is visible at the highest optical 
depths. Finally, we investigate a modification of the pair-creation threshold due to Lorentz invariance 
violation. A mild excess prevents us from ruling out an effect at the Planck energy, and we constrain 
for the first time the energy scale of the modification to values larger than sixty percent of the Planck 
energy. 

Subject headings: astroparticle physics, cosmology: observations, diffuse radiation, galaxies: active, 
gamma rays: galaxies 


1. INTRODUGTION 

The universe is not as dark as we sometimes imagine; 
even its largest voids are filled with light. The most in¬ 
tense of these photon fields, the cosmic microwave back¬ 
ground (CMB), covers the millimeter wavelength range 
and carries the relic radiation that escaped the epoch 
of recombination, less than half a million years after the 
Big Bang. At lower wavelength, from 0.1 to 1000 /rm, the 
universe is populated by the light that stars and galax¬ 
ies have emitted since the epoch of reionization {z < 10). 
Part of the light initially radiated in the ultraviolet (UV) 
and optical (O) bands is directly observable in the cos¬ 
mic optical background (GOB, 0.1 —8^m). The rest was 
absorbed by dust in the interstellar medium and around 
active galactic nuclei (AGN) and was subsequently rera¬ 
diated at lower energies, in the infrared (IR), forming the 
cosmic infrared background (GIB, 8—1000 /rm). The sum 
of the COB and GIB, the extragalactic background light 
(EBL), thus carries the 13 billion years’ radiation history 
of the universe and is a critical observable for models of 
reionization, galaxy formation and evolution, as well as 
high-energy-astrophysics phenomena, as we discuss. 

The main constraints on the EBL from observations in 
the UV-O-IR come in two flavors: direct observations, 
which tend to be contaminated by bright foregrounds 
such as the zodiacal light, and estimates from integrated 
galaxy counts, which sum t he light emitted by known 
populations of sources (e.g., Madau & Pozzetti 2000). 
The latter do not include contributions from truly dif- 
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fuse components or unobserved populations of sources, 
such as primordial stars (Pop. HI) and miniquasars that 
could have initiated the reionization of the universe (e.g. 
Madau et al. 2004 Gooray fc Yoshida|2004 ), or intra-halo 


light that was recently invoked t o explain the near- IR 
anisotropies observed by GIBER ( Zemcov et al.|[20l4 |. 

Stringent constraints also come from observations of 
gamma rays, which are more energetic than the EBL 
photons by twelve orders of magn i tude. The underly¬ 


ing proce s s, describ ed by Nikishov (1962) and Gould & 
Schreder (1967a|b), is the creation ot electron-positron 


pairs in the interaction of gamma rays from extragalactic 
sources with the EBL photon field. The survival proba¬ 
bility of a gamma ray, or gamma-ray absorption, is char¬ 
acterized by an exponential attenuation law, exp(—r), 
where the optical depth, r, depends on the redshift of 
the source and on the gamma-ray energy. 

The detection of the first distant gamma-ray sources 
led to the fi rst observational co nstraints on gamma-ray 
absorption (Stecker et al. 1992). Extragalactic sources 
observed in gamma rays, e.g. by Fermi-LAT in the high- 
energy band (HE, 0.1 — 300 GeV) or by e.g. H.E.S.S., 
MAGIC and VERITAS in the very-high-energy band 
(VHE, 0.1 — 30 TeV) are mostly AGN belonging to the 
class of blazars. The non-thermal relativistic jets emitted 
by blazars are pointed along the line-of-sight, resulting 
in an enhancement of the gamma-ray flux and energy as 
observed on Earth. Subclasses of blazars include flat- 
spectrum radio quasars (FSRQs), low-frequency-peaked 
BL Lac objects (LBLs), intermediate-frequency-peaked 
BL Lac objects (IBLs), and high-frequency-peaked BL 
Lac objects (HBLs). The energy of the peak gamma-ray 
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emission increases from FSRQs to HBLs and the UV- 
0-IR photon fields encountered by gamma rays in the 
jets appea rs to be more and mo re suppressed along this 
sequence (Ghisellini et al. 19981, making BL Lacs ideal 
cosmological probes ot the BBL (the case of FSRQs re¬ 
mains debated, see e.g. ReimerpOOT ). 

The main difficulty m constraining the EBL from 
gamma-ray observations lies in the uncertainties on the 
spectrum emitted at the source as would be observed 
on Earth without absorption, the intrinsic spectrum. If 
some curvature is present in the observed spectrum, one 
could wonder whether it is related to the emission pro¬ 
cesses occurring in the blazar’s jet or to absorption by the 
EBL during the propagation of the gamma rays to the 
Earth. Upper-limits on the EBL have been obtained as¬ 
suming no intrinsic curvature, either in the VHE or HE- 
V HE spectra of the sources , as e.g . in Meyer et al. (20121 


or 


Georganopoulos et al. (2010), and constraints naro 


bee n placed with simila r hypotheses by Orr et al. (20111 
and Sinha et al. (2014). The intrinsic HE- V H E curva¬ 
ture IS accounted for on a statistical basis in ISanchezI 


et al. (2013), through a broad-band modeling of the 
source emission from X rays to VHE in IDommguez et ah' 


( 2013 ) following a method suggested by [Mankuzhi}^ 


et al. (2010), by ex trapolating the un absorbed part of 


the HE spe ctrum in Ackermann et al. (Fermi-LAT Col¬ 
laboration, 2012), ar id by leaving the VHE curvatu re as 


a free parameter infH.E.S.S. Collaboration (2013f). In 
particular, the FerTui-LAl' and H.E.S.b. collaborations 
detected the imprint of the absorption by the EBL at 
the 6 standard deviation (a) level and at the 9tT level, 
respectively. Their approach is based on the scaling of 
existing EBL models via a normalization factor that, if 
significantly different from zero, indicates the best level of 
absorption by the EBL compatible with the gamma-ray 
data. The present work aims to overcome this model- 
dependent approach and to provide a broad-band EBL 
spectrum. 

We first show in Sec.j^that the optical depth calculated 
from the EBL spectrum can largely be simplified. One of 
the three integrals in the relation is fully reducible in an 
analytical way. We further show that the EBL spectrum 
can be deconvoluted from the pair-creation smearing ef¬ 
fect, under the assumption of a simplified evolution of 
the EBL (which we validate in Appendix [A|). Section]^ 
presents the data set s tudied. The results are discussed in 
Sec. 1^ with Sec. |4.1| focusing on the EBL spectrum and 
the room left for reionization sources or truly diffuse com¬ 
ponents. We propose a model-inde pend ent method to 
constrain the Hubble constant in Sec. 4.2 We investi gate 
the redshifts of under-constrained sources in Sec. (4.31 
And finally, we search for anomalies that could originate 
from axion-lik e par ticl es an d from Lorentz invariance vi¬ 
olation in Sec. |4.4| and |4.5[ More details on the system¬ 
atic uncertainties, best-ht parameters, and gamma-ray 
residuals are provided in appendices. 


2. EBL OPTICAL DEPTH 

The EBL optical depth to gamma rays of energy Eq 
(in the lab frame) emitted by a source at redshift zq is 


given by: 


/» 1 ^ 

J dfi ( 1 ) 


where dnjde is the density of EBL photons per infinites¬ 
imal energy band at energy e and redshift z. The other 
terms are defined in the following. 

Assuming a flat ACDM cosmology, with Hubble con¬ 
stant iLoi matter density Dm and dark energy density 
Da, the distance element is: 

dL c 1 1 _ c 

dz E[q 1 + z y^Dv + Dm(1 + z)^ Hq dz 


The integration over the energy of the EBL photons 
in Eq. is performed ab ove the pair-creation thres hold 
energy (see e.g. Eq. 8 in|Gould & Schreder||1967a|) and 
this information is encoded in the Bethe-Heitler cross 
section, with /3 S [0,1]: 


(Jryy{/3) 


3(7 J' 

TfT 


(1-/3^) 


-4/? + 2j3^ -k (3 - /3^) In 


with ctt the Thomson cross section and where 


l + P ' 
1-/3. 


( 3 ) 


2m^c^ 1 1 

Eoe 1 + z 1 — ^ 


( 4 ) 


with TOe mass of the electron and where /i is the cosine of 
the angle between the gamma ray and the EBL photons. 

Calling Co = e/(l + z) the energy of the EBL photons 
as measured today, > 0 is equivalent to: 


il + zf 


Epep 

771^ c4 


> 1 


( 5 ) 


which is the pair-creation threshold condition. 

The definition of the EBL optical depth in Eq. re¬ 
quires an integration over the distance, the energy of the 
EBL photons, and the angle between the EBL photons 
and the gamma ray. We show here that one integration 
can be reduced in a fully analytical way without any loss 
of generality. Using the change of variable^ /r —>■ /3, the 
optical depth indeed reads: 


Sarc 91 f . dn 


(1 + z)2 V Epe 


2 4 \ 2 

mic 


-P(/3max) 


( 6 ) 


^2^4 2^ 

with = 1- —1 -’ where the particle- 

EqC 1 z 

physics kernel P admits, after integration over /3, the 


^ This differs from the change of variable of|Nikishov| ||1962|| and 

.r——, . (1- '' 


Gould & Schreder 


11967a I who use /i —>■ s = (1 — 0“^) 
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following analytical expression: 


P(x) = In^ 2- ;r + 2 Li 2 

6 


1 — X 


1 — x^ 


+ (ln(l + a;) — 2 In 2) ln(l — a;) 

+ l (ln^(l-a:)-ln2(l+a:)) 


1 +3:-^ 


In 


1 


2(1-a;2) 1-a; 


(7) 


where Li 2 (a;) is the dilogarithm. 

This important result already eases the computation 
of the optical depth for a given cosmology and EBL evo- 
lution. Going further req uires an approximation: as in 
Madau & Phinney (1996), we assume that the evolution 
and spectrum ot the EBL can be locally decoupled, i.e. 

, dn, . , dn , 

de — (e, z) = deo ^ (eo , 0) X euo/ (z) (8) 

(J€ 0€q 

where evol( z) pa rametrizes the EBL history. We show 
in Appendix |A.1| that this approximation mildly impacts 
the optical depth up to redshift^of z ^ 0.6. 

Changing the variable eg —>■ Cq = ln(Aoeo/mgC‘^), the 
EBL optical depth can be written as the convolution 
product: 


t{Eq,zo) = 




En 


deo lyli, Co - In 


E. 


0 




En 


meC^ 


X lyL, (X) AT, 


X K^^{eo) 

Eq 


In 


TReC^ 


(9) 


where izli, = c j4 tt x endn(deo is the EBL specific intensity 
at z = 0 as a function of ln(eo/meC^) = ln(hz//meC^), and 
where the EBL kernel is: 


dl 

Kzoieo) =exp(-3eo) J dz — (z) 


evol{z) 

(1+^ 


X P 




Eollowing Raue fc Mazin| (2008), we parametrize the 
evolution of the EBL as evol{z) = (1 + 2;)3-/evoi^ where 
/evoi quantihes the impact of the sources. For example, 
/evoi = 0 would correspond to a simple cosmological ex¬ 
pansion of the photon held with no source term. Assum¬ 
ing a hat ACDM cosmology with Ha = 0.7, Hm_= 0.3, 
EIo = 70 km s“^ Mpc“^, we show in Appendix A.l that 


/evoi = 1. 7 best reproduces the evo lutio ns of the EBL 
model s of Franceschini et al. (2008) and Gilmore et al. 
(2012). These two state-of-the-art models from indepen¬ 
dent groups have been chosen to calibrate o ur approach, 
and we note that using a diherent model (e.g. Dominguez 
et al.||2011b| closer to /evoi = 1.2) has a minimal impact 
on our results, well below the uncertainties quoted in the 
following. 


2 This parametrization remains acceptable up to z 
optical depths smaller than 2 — 3. 


0.8 for 



Fig. 1.— EBL kernel, which yields the gamma-ray optical depth 
after convolution with the EBL intensity, as a function of the prod¬ 
uct of gamma-ray and EBL-photon energies in electron-mass units, 
in the lab frame. 


The EBL kernel from Eq. is shown in Fig. as 
a function of the squared momentum of the interact¬ 
ing photons (EBL and gamma ray), normalized to twice 
the squared electron mass. Below one, the pair creation 
threshold condition in Eq.j^is not satisfied, resulting in a 
null kernel. The kernel peaks between 2 and 4 times the 
threshold with a peak position and amplitude increasing 
with redshift. Note that for a hxed gamma-ray energy, 
the full width at half maximum of the EBL kernel almost 
spans an order of magnitude in EBL-photon energy (or 
equivalently in wavelength), whatever the redshift. This 
broadness argues against simple delta-function approx¬ 
imations to reconstruct the EBL spectrum, such as as¬ 
sumed in Sinha et al. (2014). This also implies that any 
spectral reconstruction ot the EBL based on gamma-ray 
observations yields correlated flux estimates across the 
probed wavelength range, though these correlations can 
be fully taken into account in further analyses, as demon¬ 
strated e.g. in Sec. |4.2[ 

Having shown that the EBL optical depth can be writ¬ 
ten as the convolution product in Eq. [9 we measure in 
the following the broad-band EBL spectrum based on 
local constraints and on gamma-ray absorption observed 
in gamma-ray spectra of blazars. 


3. DATASET AND ANALYSIS 
3.1. Local constraints 

Direct constraints on the intensity of the EBL from 
UV-O-IR observations come in two flavors: lower limits 
derived from galaxy counts, where faint emitters or truly 
diffuse components can be missed; and upper limits de¬ 
rived from direct measurements, which are contaminated 
by bright foregrounds (zodiacal and Galactic light), at 
least below 100/rm. In the f ollowing, we exploit the ex - 
tensive bibliographic study ofjDwek & Krennrich (|2013 ). 
We select state-of-the-art constraints from independeiit 
datasets, resulting in 27 upper limits and 25 lower limits 
(called local constraints in the following). For the latter, 
we select the estimates corrected for the completeness 
of the sample. These data are reproduced in Table 
with uncertainties combining statistical and systematic 
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Fig. 2. — Energy ranges spanned by the spectra from Tableas 
a function of redshift. 


uncertainties . In ad dition to constraints listed by |Dwek 


& Krennrich (20131, we also include the lo wer limits at 
3.6 /rm and 4.5 /rm from Ashby et al. (20131. 


procedure was followed using the 2FGL spectra (Nolan 
et al.|[20T2) for objects showing no significant variability 


in both gamma-ray bands, as long as the ground-based 
observations succeeded the launch of the satellite. 

The quasi-contemporaneity of HE and VHE data mo¬ 
tivates the use of the HE spectral shape to constrain the 
intrinsic VHE model. Photon indices have proven quite 
stable in the HE band except during flaring events, de¬ 
spite the rather high flux variability at al l time scales 


that is characteristic of blazars (see, e.g., Abdo et al. 


20101 . 


3.3. Anaysis method 
3.3.1. Test statistic 

We aim at finding the best EBL spectrum jointly ac¬ 
counting for local constraints and gamma-ray observa¬ 
tions. We define the associated test statistic as: 

= XeBL + ^ (xyray points + XhE-VHe) (H) 

7ray spectra 

X^ray points a measure of the quality of the fit of each 
gamma-ray spectrum: 


3.2. Gamma-ray data 

We study the gamma-ray spectra of blazars published 
by ground based instruments up to 2014. We have suc¬ 
ceeded in obtaining the spectral points for 106 non- 
redundant spectra, listed in Table from 38 sources. 
Representing more than 80% of the spectral data from 
blazars indexed in TeVCatj^ this is the most complete 
compilation of VHE gamma-ray data ever used for a 
study of the EBL. We paid particular attention to redun¬ 
dant data published in multiple articles, and selected the 
most up-to-date ones when a reanalysis was performed 
or when new data were included in the study. 

To perform cosmological studies, we therefore select 
the spectra of known-redshift BL Lac objects, called the 
gamma-ray cosmology sample in the following. With a 
total of ~ 270, 000 gamma rays from 86 spectra, this 
sample carries most of the information on gamma-ray 
absorption. The remaining 20 spectra, contributing an 
add itio nal ^ 30, 000 gamma rays, are used in Sections 
14.31 and l44l 

The energy ranges spanned by the gamma-ray spec¬ 
tra are shown in Fig. together with the iso-optical- 
dept h lin es derived from our best-fit EBL spectrum (see 
Sec. 4.1). The full sample covers a redshift range up 
z ~ 0.604 (PKS 1424-1-240) while the gamma-ray cos¬ 
mology sample is limited to z < 0.287 (lES 0414-1-009). 

Six oWects have uncertain redshifts, shown as red lines 
in Fig. and indicated by a question mark in the third 
colu mn of Table They are discussed in more detail in 
Sec. 14.31 

We also associated a quasi-contemporaneous HE 
gamma-ray spectrum to each VHE spectrum whenever 
available. In particular, VHE observations performed af¬ 
ter mid-2008 are contemporaneous with the sky survey 
of Fermi-LAT. Whenever HE best-fit spectral parame¬ 
ters were published in the same article as the VHE spec¬ 
trum, we stored them for further analysis. The same 


^ http://tevcat.uchicago.edu/ 


X7ray points 


E 

i G7ray points 


^ 0model(E^i; ^o) 

V 


( 12 ) 


where 4>i and are the measured flux and ass ociated 
uncertainty, and where (/^modei is discussed in Sec. |3.3.2| 
Xhe-vhe accounts for the constraints on the intrinsic 
spectral parameters. We impose, for contemporaneous 
observations in the HE and VHE bands, that the intrinsic 
VHE spectrum be softer than the HE spectrum: 


Xhe-vhe — ©(rHE — B) 


The — T 

^rHE 


(13) 


where Fhe and oThb are the photon index and associated 
uncertainty measured at HE, and where F is the intrinsic 
VHE photon index. For curved spectra, the indices and 
uncertainties are computed at the intersection of the HE 
and VHE range. 

Xebl accounts for the lower limits and upper limits 

listed in Table 1 1 Only the constraining 

sides of the limits are considered, using again Heavyside 
functions: 

Xebl = E ~ ^Bj/(Ai)) x 

i G LL 

+ y:e(.c(A.)-Ox(AbJzjA)' 

iGUL A u J 

(14) 

where vA(A) is the model of the data, and where LL 
and UL denote the ensembles of lower and upper limits, 
respectively. 

The local constraints effectively taken into account in 
Eq. depend on the EBL-wavelength range probed by 

the gamma-ray data. The threshold condition in Eq. 
imposes A = /ic/eo < Amax, where Amax is the maximum 
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TABLE 1 

Local constraints on the EBL spectrum used in this p aper, largely extracted from|Dwek & | 

Krennrich| (|2013[). 


A 

Lower limit 

Upper limit 

Experiment 

Reference 

/im 

nW m~^ sr~^ 

nW m~^ sr“l 




0.153 

1.03 ±0.15 


Galex 



jXu et al. 


2005 

0.1595 

3.75 ± 1.25 


HST/WFPC2 

G 

ardner et al. 


^UUU 

0.231 

2.25 ±0.32 


Galex 



Au et al. 


yUUb 

0.2365 

3.6 ±0.5 


HST/WFPC2 

G 

ardner et al. 


2000 

0.3 


18 ± 12 

HST/WFPC2 



iBernstem 


^UUV 

0.3 

3.7 ± 0.7 


HST±ground 


1 

Tbtani et al. 


2001 

0.4 


26 ± 10 

dark cloud 




Mattila 


lyyo 

0.44 


7.9 ±4.0 

Pioneer 10/11 

iMatsuo 

<a et al. 


^Uil 

0.45 

6.1 ± 1.8 


HST±ground 


1 

iotani et al. 


^UUl 

0.5115 


30 ±9 

ground 



liJube et al. 


lyyy 

0.55 


55 ±27 

HST/WFPC2 



Ikiernstein 


^UUV 

0.61 

7.4 ± 1.5 


HST±ground 


J 

Tbtani et al. 


2001 

0.64 


7.7 ± 5.8 

Pioneer 10/11 

liVlatsuolca et al. 


2011 

0.81 

9.3 ± 1.6 


HST±ground 



iotani et al. 


2DOT 

0.814 


57± 32 

HST/WFPC2 



lisernstem 


yuuv 

1.25 


21 ± 15 

COBE/DIRBE 

iLevenson et al. 


2007 

1.25 

11.5 ± 1.3 


HST±ground 


3 

iotani et al. 


SDUT 

1.25 

11.7± 2.6 


Subaru 


ixeenan et al. 


l^UiU 

1.6 

11.5± 1.5 


Subaru 


Keenan et al. 


^UiU 

2.12 

10.0 ±0.8 


Subaru 


ixeenan et al. 


smn 

2.2 


20 ±6 

COBE/DIRBE 

1^ 

evenson et al. 


2007 

2.2 

9.0 ± 1.2 


HST±ground 



iotani et al. 


^UUl 

3.5 


13.3 ± 2.8 

COBE/DIRBE 

1^ 

e-' 

/enson et al. 


2DD7 

3.6 

5.6 ± 1.0 


Spitzer/IRAC 



Ashbv et al. 


2013 

4.5 

4.4 ±0.8 


Spitzer/IRAC 



Ashbv et al. 


2013 

4.9 


22 ± 12 

COBE/DIRBE 

lAre 

ndt 66 Uwek 


2003 

15 

1.9 ±0.4 


ISO/ISOCAM 

ItioDwooa et al. 



16 

2.2 ±0.2 


Spitzer 


l iepiitz et al. 


^Uil 

24 

2.86t°;l® 


Spitzer/MIPS 

Bethermin et al. 


2010 

60 


28.1 ± 7.2 

COBE/DIRBE 

f mkbeiner et al. 


^UUU 

65 


12.5 ±9.3 

Akari 

iMatsuura et al. 


^Uil 

70 


Spitzer/MIPS 

Bethermin et al. 


2010 

90 


22.3 ±5.0 

Akari 

IMatsuura et al. 



100 

8.35 ±0.95 


Herschel/PACS 



liJerta et al. 


2011 

140 


12.6 ±6.0 

COBE/FIRAS 


K ixsen et al. 


lyyy 

140 


20.1 ± 3.6 

Akari 

Iiviatsuura et al. 


2011 

140 


15.0 ±5.9 

COBE/DIRBE 

lUdegara et ai. 


^UUV 

160 


Spitzer/MIPS 

Bethermin et al. 


2010 

160 


13.7 ± 6.1 

COBE/FIRAS 

It' ixsen et al. 


■rags 

160 


13.7 ± 4.0 

Akari 

IMatsuura et al. 


2011 

160 


14.4 ± 2.4 

Spitzer/MIPS 



i'enin et al. 


'2iJV2 

240 


10.9 ±4.3 

COBE/FIRAS 


J 

t ixsen et al. 


■rags 

240 


12.7 ± 1.6 

COBE/DIRBE 

ludeeara et ai. 


2007 

250 

10.13l|®o 


Herschel/SPIRE 

Bethermin et al. 


2012 

250 


10.3 ±4.0 

COBE/FIRAS 


li^ ixsen et al. 


lyyy 

350 

6A6tl:lt 


Herschel/SPIRE 

Bethermin et al. 


2012 

350 

2.80lg;®? 

5.6 ± 2.1 

COBE/FIRAS 

It' ixsen et al. 


■rags 

500 


Herschel/SPIRE 

Bethermin et al. 


2012 

500 


2.4 ± 0.9 

COBE/FIRAS 


K ixsen et al. 


lyyy 

850 

0.24 ±0.03 


SCUBA 



emcov et al. 


l^UiU 

850 


0.50 ±0.21 

COBE/FIRAS 


li^ ixsen et al. 


iyy8 


wavelength associated with each spectrum, with: 

x{l + zof (15) 


, _ ^ ^'max.VHE 

^max — ^ o 

rrieC 


95 X 


rrieC 
^^max,VHE 

20TeV 


X (1 + 


imization over the parameters of the intrinsic spectra. 
This minimization is performed for each spectrum using 
the MIGRAD method of MINUIT. The minimization over 
the EBL parameters is performed in three steps, using 
successively SIMPLEX, MIGRAD, and HESSIAN. 

3.3.2. Intrinsic spectral models 


The minimum wavelength is set a posteriori by succes¬ 
sively adding free parameters to the model until no fur¬ 
ther improvement is found, as discussed in Sec. 3.3.3| 


The VHE gamma-ray spectra are modeled with 

4>mode\{E, z) = 4>iryt(.E) X exp {-t{E, z)) (16) 


2013f|, the minimization of the in Eq. 11 

in turn as done by the 

iUUCi, iO ClCLCiiii 

H.E.S.S. Collaboration 

(2013f 


As a 
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sen: a power law (PWL in Table HI). 4)'p^n'l{E) = 
(I)q{E/Eq)~^, where Eq is the reference energy, 
is the flux normalization and P is the photon in¬ 
dex. The reference energy is fixed, in this work, to 
the central value of the energy range of each spec¬ 
trum, i.e. Eo = yj i^min,VHE-Emax,VHE• To Search for 
intrinsic curvature, we also test the three-parameter 
log parabola (LP), (l>i,v{E) = 

where a is the photon index at Eq and where b is 
the curvature parameter, the exponential cutoff power 
law (EPWL), 0 epwl(£)) = ME/Eq)-'^ exp{-E/Ecnt), 
where E^ut is the cutoff energy. We also consider 
an exponential cutoff log parabola (ELP), (/)elp(E) = 
exp(—E/Ecut)) however none of 


the spectra is found to signihcantly prefer this model. We 
impose concave intrinsic spectra, i.e. b > 0 and Ecut > 0. 

For a given set of intrinsic models, we inspect the resid¬ 
uals of the individual spectra fixing the EBL parame¬ 
ters to their best-fit values. If a more complex intrinsic 
model is preferred at least at the 2a level, it is chosen for 
the next iteration, until the set of models converges. To 
avoid any overrestriction of the parameter space, once 
the stable set of intrinsic models was found, we suc¬ 
cessively changed each intrinsic spectrum to its more 
complex counterpart, whichever was most preferred, and 
checked that it had no impact on other models nor on 
the best-fit EBL spectrum. This ensured the robustness 
of our method and of the choice of the intrinsic models. 


TABLE 2 

Gamma-ray spectra used in this paper. 


Source 

Class 

Redshift 

Experiment 

Obs. Period 

Model 




Reference 



IC 310 

HBL 

0.019 

MAGIC 

2009-2010 (low) 

PWL 



Aleksic et al. 

2014b1 




MAGIC 

2009-2010 (high) 

PWL 



Aieksic et al. 

2U14bi 

Markarian 421 

HBL 

0.031 

CAT 

1998 

PWL 

IPiron i! 

tJAi uoiiab 

■ 

2DDDI 




HEGRA 

1999-2000 

PWL 


Anaronian et ai 


2UU2I 




HEGRA 

2000-2001 

LP 


Anaronian et ai 


SDTOi 




I'ibet 

2000-2001 

PWL 


Amenomori et al 

' 

2UU3I 




HESS 

2004 

EPWL 


Aharonian et al. 


OUbci 




MAGIC 

2004-2005 

EPWL 



Ij 

Albert et al. 


OOYd 1 




MAGIC 

2006-04-22 

PWL 




Aleksic et al 

■ 

2010 1 




MAGIC 

2006-04-24 

PWL 




Aleksic et al 


20101 




MAGIC 

2006-04-25 

PWL 




Aleksic et al 


20101 




MAGIC 

2006-04-26 

PWL 




Aleksic et al 


20101 




MAGIC 

2006-04-27 

EPWL 




Aleksic et al 


20101 




MAGIC 

2006-04-28 

PWL 




Aleksic et al 


20101 




MAGIC 

2006-04-29 

PWL 




Aleksic et al 


20101 




ARGO-YBJ 

2007-2010 (fluxl) 

PWL 




Bartoli et al 


YlUii 1 




ARGO-YBJ 

2007-2010 (flux2) 

LP 




Eartoli et al 


20111 




ARGO-YBJ 

2007-2010 (flux3) 

EPWL 




Bartoli et al 


YlUii 1 




ARGO-YBJ 

2007-2010 (flux4) 

PWL 




Bartoli et al 


YlUii 1 




VERITAS 

2008 (very low) 

LP 




Lcciari et al. 

2 

011b 1 




VERITAS 

2008 (low) 

LP 



Acciari et al. 

2Ullbi 




VERITAS 

2008 (mid) 

LP 



Acciari et al. 

2Ullbi 




VERITAS 

2008 (high A) 

LP 



Acciari et al. 

2Ullbi 




VERITAS 

2008 (high B) 

EPWL 



Acciari et al. 

2Ullbi 




VERITAS 

2008 (high C) 

PWL 



Acciari et al. 

2011b 1 




VERITAS 

2008 (very high) 

EPWL 



Acciari et al. 

2Ullbi 




TACTIC 

2005-2006 

PWL 



i^iinarma et al 

■ 





TACTIC 

2008 

PWL 



unanara et ai 


20101 




TACTIC 

2009-2010 

PWL 



Chandra et al 


’zurz ) 

Markarian 501 

HBL 

0.034 

HEGRA 

1997 

PWL 


lAn 

aronian et al 


20011 




TACTIC 

2005-2006 

PWL 


iLroaambe et ai 


y:uu» I 




MAGIC 

2006 

PWL 

lAndernuo et ai. 

‘J 

OOObi 




ARGO-YBJ 

2008-2011 

PWL 




Eartoli et al 

■ 

20121 




ARGO-YBJ 

2011 (high) 

PWL 




Bartoli et al 


’^urz ) 

lES 2344-1-514 

HBL 

0.044 

Whipple 

1995 (dataset B) 

PWL 


jSch 

roedter et al 

' 

2005 1 




MAGIC 

2005-2006 

PWL 



{Albert et al. 


007c 1 




VERITAS 

2007 (low) 

LP 



Acciari et al. 

2Ullai 




VERITAS 

2007 (high) 

PWL 



Acciari et al. 

2Ullai 




MAGIC 

2008 

PWL 



lAleksic et a 

■ 

20131 

Markarian 180 

HBL 

0.045 

MAGIC 

2006 

PWL 




{Albert et al 

' 

2DD0i 

lES 1959-1-650 

HBL 

0.048 

HEGRA 

2002 

PWL 

Ah 

aronian et al. 


003a 1 




Whipple 

2002 

PWL 




luaniei et a 


YIUUO 1 




MAGIC 

2006 

PWL 


1 laeiiaterri et ai 


2DDSi 




VERITAS 

2007-2011 

PWL 




lAiiu et al.l 

‘J 

013b 1 

BL Lacertae 

IBL 

0.069 

MAGIC 

2005 

PWL 



lAibert et al.l 

2DD7ai 




VERITAS 

2011 

PWL 




lArien et ai 

T|M3i 

PKS 2005-489 

HBL 

0.071 

HESS 

2004-2007 

PWL 

|H.E.S.S. Collaboration! 

2UlUb 1 

RGB J0152+017 

HBL 

0.08 

HESS 

2007 

PWL 

lAharoman et al.l 

2U08a 1 

SHBL J001355.9-185406 

HBL 

0.095 

HESS 

2008-2011 

PWL 

jH.E.S.b. uoiiaoorationl 

2U13ai 

W Comae 

IBL 

0.102 

VERITAS 

2008-01-04 

PWL 



lAcciari et al 

iISddsi 

lES 1312-423 

HBL 

0.105 

HESS 

2004-2010 

PWL 

H.E.S.b. uoiiaborationj 

2013d 1 

VER J0521-I-211 

IBL 

0.108 

VERITAS 

2009-2010 

PWL 

lArcUambault et al 

■ 

20131 

PKS 2155-304 

HBL 

0.116 

HESS 

2002-06 

PWL 


Aharonian et al. 

2 

OObai 




HESS 

2002-10 

PWL 


Aharonian et al.l 

wnuai 




HESS 

2003-06 

PWL 


AEaroman et al.| 

2UUba 1 
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TABLE 2 — Continued 


Source Class Redshift Experiment Obs. Period Model Reference 





HESS 

HESS 

HESS 

HESS 

HESS 

HESS 

MAGIG 




HESS 

HESS 

HESS 

B3 2247-1-381 

HBL 

0.119 

MAGIC 

RGB J0710+591 

HBL 

0.125 

VERITAS 

H 1426-1-428 

HBL 

0.129 

HEGRA 

HEGRA 

lES 0806-1-524 

HBL 

0.138 

VERITAS 

lES 0229-1-200 

HBL 

0.14 

HESS 

VERITAS 

IRXS J101015.9-311909 

HBL 

0.143 

HESS 

H 2356-309 

HBL 

0.165 

HESS 

HESS 

HESS 

RX J0648.7+1516 

HBL 

0.179 

VERITAS 

lES 1218+304 

HBL 

0.182 

VERITAS 

VERITAS 

lES 1101-232 

HBL 

0.186 

HESS 

lES 0347-121 

HBL 

0.188 

HESS 

RBS 0413 

HBL 

0.19 

VERITAS 

lES 1011+496 

HBL 

0.212 

MAGIC 

lES 1215+303 

HBL 

0.237? 

VERITAS 




MAGIC 

S5 0716+714 

IBL 

0.26? 

MAGIC 

PKS 0301-243 

HBL 

0.266 

HESS 

lES 0414+009 

HBL 

0.287 

HESS 

VERITAS 

3C 66A 

IBL 

0.335? 

MAGIC 

PKS 0447-439 

HBL 

0.343? 

HESS 

PKS 1510-089 

FSRQ 

0.361 

HESS 

4C +21.35 

FSRQ 

0.432 

MAGIC 

PG 1553+113 

HBL 

0.433? 

HESS 

MAGIC 

MAGIC 

MAGIC 




VERITAS 




HESS 

3C 279 

FSRQ 

0.536 

MAGIC 

PKS 1424+240 

IBL 

0.604? 

VERITAS 

VERITAS 

VERITAS 




MAGIC 

MAGIC 

MAGIC 


2003-07 

PWL 


Aharonian et al. 

( 

2005a 

2003-08 

PWL 


Aharonian et al. 


20Uba 

2003-09 

PWL 


Aharonian et al. 

■ 

20Uba 

2003-10-11 

PWL 


Anaronian et ai. 

1 

2UUbb 

2005-2007 

PWL 

H.E.S 

.S. 

UoiiaDoration 


ZUiUC 

2006-07/08 

LP 


Abramowski et a] 


■ 


2006-07/08 

LP 



Aieksic et al. 



0I2c 

2006-08-02 

PWL 

H.E.S.S. 

>^oiiaDoration 


zuiza 

2006-08-03 

PWL 

M.ii.b.Si. uoliaDoration 


:^ui:^a 

2008-08-09 

PWL 


lAUaroman et a 


(| 2 uuy 

2010 

PWL 



Aieksic et al. 


2012a 

2008-2009 

PWL 



Acciari et al. 

1 

20iUb 

1999-2000 

PWL 


Ah 

aronian et al. 

1 

2UU3b 

2002 

PWL 


Aharonian et al. 


20U3b 

2006-2008 

PWL 



Acciari et al. 


20Uya 

2005-2006 

PWL 


[Ah 

aronian et al. 


2UUVc 

2009-2012 

PWL 




lAiiu et a 


(|2UI4 

2006-2010 

PWL 

H.E.S.S. 

UoiiaDoration 


2012c 

2004 

PWL 

n.ii.b.b. uoiiaDoration 


zuiua 

2005 

PWL 

n.ii.b.b. UoiiaDoration 


2UIUa 

2006 

PWL 

H.B.S.S. Collaboration 


:^uiua 

2010 

PWL 




lAiiu et a 


(|2UIi 

2007 

PWL 



Acciari et al. 


2009 b 

2008-2009 

PWL 



Acciari et al. 


2010a 

2004-2005 

PWL 


Ah 

aronian et al. 


20UVa 

2006-09-12 

PWL 


Anaro 

man et al. 


20UVb 

2008-2009 

PWL 




Aliu et al. 


2012a 

2007 

PWL 



|A 

bert et al. 


200Vb 

2008-2012 

PWL 



Iaiiu et al. 


2013a 

2011 

PWL 


lAieKsic et al. 


2012b 

2007-2008 

PWL 


AnaernuD et ai. 


20Uya 

2009-2010 

PWL 

H.E.S 

.S. UoiiaDoration 


2U13c 

2005-2009 

PWL 

ti.ii.Jb.Si. UoiiaDoration 

1 

20i2b 

2008-2011 

PWL 



lAiiu et ai. 

1 

20i2b 

2009-2010 

PWL 



Aleksic et al. 

1 

20IIb 

2009-2010 

PWL 

H.E.S.S. 

.;oiiaDoration 

1 

zuijb 

2009 

PWL 

H.E.S.S. 

Collaboration 


zuide 

2010-06-17 

PWL 


lAleicsic et al. 


20IIa 

2005-2006 

PWL 


Ah 

aronian et al. 


2008b 

2007 

PWL 



Aleksic et al. 

1 

20I2d 

2008 

PWL 



Aleksic et al. 

1 

20I2d 

2009 

PWL 



Aleksic et al. 

1 

2 

mm 

2010-2012 

PWL 




lAliu et a] 


■ 

2015 

2012 

PWL 


Abramdwski et al. 


2015 

2006 

PWL 

IMA 

jiu UoiiaDoration 


2DDS 

2009 

PWL 

'-7 

Vrchambault et al. 


2014 

2011 

PWL 

7 

Vrchambault et al. 


2014 

2013 

PWL 

7 


lamDauit et a 


' 

mm. 

2009 

PWL 



Aleksic et al. 


'2 

0I4a 

2010 

PWL 



Aleksic et al. 


2014a 

2011 

PWL 



Aleksic et al. 


2014a 


3.3.3. Deconvolution method 

Similarly to the splines used in|Mazin & Raue| (|2007|), 
we chose a spectral model for the EHL that is a lin- 
ear function of its parameters. As shown in Eq. the 
pertinent variable is the logarithm of the EBL-jmoton 
energy, or equivalently the logarithm of the EBL wave¬ 
length I = ln(A/Aref), where Aref is a reference wave¬ 
length, set e.g. to 1 fiui. We model the EBL intensity 
as a sum of Gaussians of this variable, with fixed widths 
and peak positions: 


T ^n f ~ 

exp I-1 

i ^ I / 

= y^aiJ\f{l;li,ai) (17) 

i 

The peak positions are logarithmically spaced: 

l0g(Ai/Aref) = li = lo+iX lS .1 (18) 

The width of the Gaussian functions depends on the 
distance between anchor points A/. We impose that the 
sum of two consecutive Gaussians of unity amplitude be 
equal to one right in between the two anchor points, i.e. 

Ai = aix 2\/2 In 2 (19) 

The Gaussians leak into the neighboring bins and the 
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average flux in a bin centered on li is then 

^4 = 2 ^ 0 ,— / dZexp- 

j Jh-Ai/2 \ 


( 20 ) 


with 


dVln^ 


erf ^(z — j + ^) X 2Vln2^ 
— erf M * — J — -) X 2Vln 2 


( 21 ) 


The parameters are left free to vary. Once best-fit 
values, uncertainties, and correlations have been deter¬ 
mined, the binned average flu x vl l can easily be derived 
from the linear relation in Eq. 20 This linearity permits 
the propagation of uncertainties for Gaussian distribu¬ 
tions of the weight Oi, fully accounting for the correlation 
ter ms. Th is justifies the use of the HESSIAN method in 
Sec. |3.3.1[ which yields such symmetric Gaussian uncer¬ 
tainties. 

Given the linearity of the model with respect to its 
parameters Ui, the optical depth can be rewritten as: 


t(Eo,zo) =Y°'^ U{Eq,zq) 


( 22 ) 


where 


, ,-1^ \ dTrCry Eq 

ti{Eo,zo) = X 


X A/'(-; 0, ai) ® ( e* -I- In 


he/Xi 


Et) 


(23) 


where e,- = In 

One can compute the weights ti(Eo,zo) in the very 
beginning of the fitting procedure, further reducing the 
computation expense. For a set of about 90 spectra and 
associated models, the full fitting procedure of the EBL 
spectrum takes about ten seconds of CPU time on a 
3 GHz core, highlighting the significance of the analyt¬ 
ical work shown in this section and in Sec. [21 

4. RESULTS 
4.1. EBL spectrum 

The best-fit spectral models of the intrinsic spectra are 
listed in column 6 of Table Most of the spectra of the 
gamma-ray cosmology sam^e (71/86) are best described 
by PWL models. The other fifteen spectra, modeled by 
LP and EPWL models, correspond either to intensive 
campaigns in low states of prominent sources (2004-05 
MAGIC campaign of Markarian 421, large-zenith-angle 
H.E.S.S. observations of Markarian 421, 2007 VERI¬ 
TAS campaign on lES 2344-1-514, 2007-2010 campaign 


of ARGO-YBG on Markarian 421: |Albert et al.| 

2007dl 

Aharonian et al.|2005cl 

Acciari et al.|2011a||Barto 

li et al. 

2011 respectively), to observations of tlares (hig 
ot Markarian 421 observed by HEGRA in 2000-20 

1 state 
01 and 


MAGIC in 2004-2005, major outburst of PK S 2155-304 
in 2006 observed by MAGIC and H.E.S.S.: Aharonian 



Fig. 3.— EBL intensity at 2: = 0 as a function of wavelength. 
The best-fit spectra derived in this work are shown with light 
blue (gamma rays only, four-point spectrum) and blue points 
(gamma rays -j- direct constraints, eight-point spectrum). Lower 
and upper limits are shown with orange upward-going and dark- 
brown down ward-going arrows , respe ctivel y. For comparison with 
the wo rk of|Ackermann et al.| (|2012|| andjH.E.S.S. Collaboration! 
( 2013f||' the 1 a (stat -j- sys) confour of the oest-tit scaled-up model 
1 Uilmore et al. 12012^ is shown as filled blue region, using a scaling 
factor ot 1.L6 as snown in Table [4| 


et al.|20d2[ [Aleksic et al. 


]|2010| |20I2c| [Abramowski ^ 

2013[ respectively), or to both (2006-2008 campaign on 
Markarian 421 by VERITAS, including flares: Abdo et al. 
20111. In such cases, enhanced statistics at the highest 
energies enable the probe of intrinsic curvature. No spec¬ 
trum is preferentially modeled with an ELP model. 

The 86 spectra probe the wavelength range 0.26 — 
105 /xm, for a bin size A; = 0.75. The maximum 
wavelength corresponds to the pair-creation threshold, 
as described in Eq. |15| A smaller minimum wave¬ 
length would result in an underconstrained EBL inten¬ 
sity in the first bin. A smaller binning does not sig¬ 
nificantly improve the quality of the fit. With a total 
of 630 points and 187 free parameters for the intrinsic 
spectra, the best-fit model results in a test statistic of 
(Ex^ray points + Xhe-vhe)Ac(/ = 340.1/443. The small 
value of the reduced is not surprising, as the corre¬ 
lations between gamma-ray spectral points are not ac¬ 
counted for when fitting such archival data (the gamma- 
ray community is only starting to publish covariance ma¬ 
trices for spectral points). The uncertainties are also as¬ 
sumed to be Gaussian (underlying assumption for the 
test), while a full treatment at the event level would 
account for the Poisson statistics of the events from back¬ 
ground and signal regions. 

The constraint from the hardness of the HE spectra as¬ 
sociated with the VHE observations, proves a posteriori 
to play a minor role, y) Xhe-vhe <0-1, indicating that 
there is no tension with the assumption of broad-band 
concavity in the intrinsic spectra. No tension is found 
with the local EBL constraints either, with Xebl — 2-4 
to which uebl = 7 local constraints contribute. Both 
the local EBL constraint and the hardness constraint 
thus barely impact the best-fit estimate of the EBL spec¬ 
trum, but they nonetheless play a significant role when 
the spectrum departs from the best-fit point, thus im- 


et al.l 



















































EBL, Hubble constant, and anomalies 


9 


pacting the uncertainties on the EBL. For the binning 
and the wavelength range probed here, the Gaussian-sum 
model admits eight parameters, resulting in a total test 
statistic of x^/nd/ = 342.5/442. The number of degrees 
of freedom accounts for eight free parameters to model 
the EBL and seven local points constraining this model, 
in addition to the 443 degrees of freedom accounting for 
the gamma-ray spectra. 

Using the same intrinsic models and a null ab¬ 
sorption, the gamma-ray spectra best-fit test statis¬ 
tic is E(X7ray points + xIe-VHe) = 489.1. _ With 
eight additional free parameters, the Gaussian-sum 
model is preferred by the gamma-ray data at the 


\/2 erfc = 489.1 — 340.1)] = Her level, where 

Vs is the x^ probability for eight degrees of freedom. 
This signihcance level may not readily be interpreted as 
the detection significance of the EBL signature, because 
the intrinsic models chosen for the likelihood ratio test 
are based on the best-fit EBL spectrum. Conservatively 
choosing instead the intrinsic models based on fits di¬ 
rectly to the observed spectra assuming no EBL absorp¬ 
tion results in eight intrinsic models being changed to a 
more complex or an equally complex model. The best- 
fit EBL model is then preferred at the 10.3 a level, with 
no significant impact on the best-fit parameters. The 
cumulative EBL effect therefore leaves a strong imprint 
in the gamma-ray spectra, at a high significance level in 
between the two values given above. The EBL intensi¬ 


ties computed from the best-ht parameters, as in Eq. 20 
are listed in Table [3 and shown in Fig. (blue points). 
The overall systematic uncertainty on tne EBL level is 
estimated in Appendix [A| to be on the order of 7 — 8 %. 


TABLE 3 

Best fit EBL intensities. 


A 

/i.m 

-^min 

Amax 

/im 

Viv 

nW sr~^ 

0.38 

0.26 

0.55 

7.4 ±3.7 

0.80 

0.55 

1.2 

12.4 ± 1.9 

1.7 

1.2 

2.5 

10.8 ±0.8 

3.6 

2.5 

5.2 

6.7 ± 0.7 

7.6 

5.2 

11 

0.95 ±0.53 

16 

11 

23 

3.79 ±0.26 

34 

23 

50 

3.14 ±0.25 

72 

50 

105 

9.8 ± 1.1 


Theoretical and empirical EBL intensities can sim¬ 
ilarly be compared to local constraints and gamma-ray 
data. Table |4] lists the test statistics obtained with a 


(2010) (kUlU). Al. 


null EBL, the best-ht eight-point EBL spectrum, and 

the models of Gilmore et al. (|2012|) (G12), Franceschini 

et al. (20081) (k08). 

Dominguez et al. (2011bD (Dll), 

ytecker et al. (2006 

(SU6) khaire & Srianand (|2015|) 

(KSi4),|Einkc! k al. 

(20101 (ElO), andjKneiske & Dole 


the EBL models are preferred at 


the 10 — 12(7 level to a null EBL (column 7), as com¬ 
puted from the improvement in the quality of the fit of 
the gamma-ray data (column 4). The EBL models can 
also be compared to the best-ht eight-point spectrum 
with a likelihood ratio test, assuming that the models 
are nested. Such an approach is justihed by a differ¬ 
ence in optical depth between a model and its eight-point 
Gaussian-sum ap proxi mation smaller than 2 % (see Ta¬ 
ble]^ of Appendix|A.2 ). The difference between the best- 


ht eight-point spectrum and the models of G12, F08, and 
Dll corresponds to a small 0.5 — 1.5 a preference for the 
Gaussian-sum description (column 8). The EBL shapes 
predicted by S06, KS14, FlO, and KIO are disfavored at 
the > 3 — 6 (7 level with respect to our be st-ht spectrum. 
We have not included here the models of|Helgason et al.| 
(|2012| and 
only up to 


Stecker et al. (2012|), which estimate the EBL 
he M (~ 5 fim) and I (^ 0.8 /rm) photometric 
bands, respectively, and are thus not relevant for most 
of the TeV spectra listed in Table We note nonethe¬ 
less that their results are within 1 — 2 cr of the estimates 
derived in the present work at small wavelengths. Follow¬ 
ing the approach of the EB L studies by the H.E.S.S. anti 
Fermi-LAT collaboratio ns (|H.E.S.S. Collaboration 2013f 
Ackermann et al. , we estimate the best-ht normal- 


ization factor of the EBL models, as listed in column 9 
The best-ht normalizations of the models listed in Table 
are all larger than 1, though only at the 1 — 2a level for 
the models of G12, F08, Dll, and S06. This is not sur¬ 
prising as most of these models were designed to match 
sets of lower limits based on less complete surveys than 
the most recent studies, and as such tend to yield lower 
estimates of the EBL intensity. We compare in Fig.[^the 
scaled-up model of G12 with constraints derived from the 
measurements of H.E.S.S. below z < 0.2 and Fermi-LAT 
between 0.5 < z < 1.6. The H.E.S.S. and Fermi-LAT 1 a 
conhdence contours include the statistical and systematic 
uncertainties on the measurement of the parameter nor¬ 
malizing the model of F08 . in the wavelength ranges rel¬ 
evant to thes e studies (seejH.E.S.S. Collaboration] 2013f 
Biteau||2013l). A good agreement on the level of EBL is 


found between these works and our results. 

We investigate more closely the origin of the differ¬ 
ences between the gamma-ray based EBL spectrum and 
the models. Good agreement between our results and all 
of the models in TableHis found below 5 /im. The models 
of S06, KS14, FlO, antl^DlO are thus disfavored because 
of the rather low level of GIB they predict. The mod¬ 
els of G12, F08, and Dll tend to predict a level of EBL 
that is higher than the gamma-ray estimate around 8 fim, 
while lower around 16 /xm. The ~ 16 fxm point from the 
gamma-ray data is found to be ~ 3(7stat+sys above the 
galaxy-count estimate, but the lack of an observation 
around 8 ^m prevents a direct comparison at this wave- 
lengtl0 These hints of deviations in the 5—20 /im region, 
which could be indicative of the signature of polycyclic 
aromatic hydrocarbons, should be taken with caution as 
neighboring points from the gamma-ray estimates are 
correlated. Using the covariance matrix shown in Ap¬ 
pendix]^ the correlation coefficient between the ~ 8 /im 
and ^ 16/im points is negative (—16%), showing that 
a decrease in the latter would result in an increase of 
the former, which would considerably reduce the devia¬ 
tions. A hrm conclusion on the signature of polycyclic 
aromatic hydrocarbons cannot be drawn at this stage, 


Renormalizing theoretical or empirical EBL mociels is not phys¬ 
ically motivated, since different components contributing to the 
spectrum need not scale in the same way. Nonetheless, measuring 
the best-fit scaling factor for a model is a straightforward way to 
quantify the model’s compatibility with the experimental gamma- 
ray data. 

® A linear interpolation in log(A) between the 4.5 and 16 ^m 
points from galaxy counts yields a level of EBL ~ 3 iTstat+ays higher 
than the gamma-ray estimate. 
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as indicated by the small difference in quality between 
the model-independent fit of this work and the model- 
dependent fits (column 2 in Table |^. Increased gamma- 


ray statistics from local sources (z ^ 0.1) and improved 
constraints from galaxy counts in the 5 — 20 fj,m region 
would definitely help in deciding the matter. 


TABLE 4 

Test statistics of the EBL models. 


Model 


ndf 

Xj 

^EBL 

liEBL 

cr{- ^ no EBL) 

(T(this work 7 ^ •) 

^ a 

X^(«) 

Aebl(“) 

liBBLfo) 

No EBL 

489.1 

443 

489.1 









This work 

342.5 

442 

340.1 

2.4 

7 

11.0 






G12 

353.5 

455 

345.2 

8.4 

12 

12.0 

0.5 

1.13 ±0.07 

346.0 

1.9 

4 

F08 

366.1 

451 

350.2 

15.9 

9 

11.8 

1.3 

1.05 ±0.07 

352.1 

13.6 

6 

Dll 

370.0 

453 

351.0 

19.0 

10 

11.8 

1.5 

1.16 ±0.05 

356.9 

3.5 

5 

S06 

392.7 

445 

382.0 

10.7 

3 

10.3 

5.0 

1.18 ±0.08 

380.0 

7.5 

3 

KS14 

401.4 

456 

362.1 

39.3 

13 

11.3 

3.0 

1.44 ± 0.07 

361.2 

0.5 

4 

FIO 

424.0 

449 

390.2 

33.8 

7 

9.9 

5.7 

1.48 ±0.07 

384.6 

3.6 

3 

KDIO 

433.0 

457 

366.7 

66.3 

14 

11.1 

3.5 

1.52 ±0.14 

346.5 

0.2 

1 


Significance with which the gamma-ray data prefer the 
model to the absence of EBL. 

Significance with which the gamma-ray data prefer the 
EBL spectrum from this work to the model. 


To estimate the total brightness of the EBL, we per¬ 
form a fit similar to that described above, extending 
the wavelength range to 0.1 — 1000 p,m and using twelve 
points, corresponding roughly to A/ ~ 0.75. The fit is 
underconstrained below 0.25/tm and the EBL estimate 
in the lowest wavelength bin is compatible with zero at 
the 1 cr level. The brightness of the COB in 0.1 — 8 /im is 
estimated to be 36 ± llnW m“^ sr“^, with a large un¬ 
certainty coming from the region A < 0.25 ^m. Above 
105 ^m, the gamma-ray data do not constrain the fit but 
the CIB is tightly constrained by local measurements, 
yielding a 8 — 1000 /rm CIB brightness compatible with 
that of the COB, with a value of 25.9±3.4 nW m“^ sr“^. 
The total brightness of the EBL is measured with a 
^ 20% accuracy as 62± 12 nW m“^ sr“^, which is equiv¬ 
alent to 6.5 ± 1.2% of the brightness of the CMB. This 
is compatible at the 1 a level with previ ous estimations 
based on galaxy counts (|Dole et al.|2006D and on models 
(see e.g. Table 1 in G12, Table 4 in Ull, and Table 2 in 
FIO). 

We also show in Fig. the best-fit EBL spectrum 
derived taking into account only the gamma-ray data 
(light-blue points). A larger binning (A; = 1.5) results 
because of the loss of information. Gamma-ray data show 
a spectrum of the EBL from mid-UV up to far IR that 
is in good agreement with estimates based on galaxy 
counts. The Her deviation from a null EBL intensity 
and the fact that EBL reconstructed from gamma-ray 
observations follows the expected spectrum are pieces of 
evidence against scenarios in which the VHE emission 
of blazars would primarily come from ultra-high-energy 
cosn iic rays (UHEGR) repro cessed along the line of sight 
(e.g. Essey fc Kusenko| 20101. Quite a cosmic conspiracy 
would be needed to explain how secondary gamma rays 
from UHEGR carry the very same imprint as that ex¬ 
pected from absorption of primary gamma rays. A more 
quantitative study of the impact of this process on the 
EBL reconstruction is left for future work. 

To quantify the compatibility between galaxy counts 
and the four-point EBL spectrum based only on gamma- 
ray data (light-blue points in Fig.|^, we compute the dif¬ 
ference between these estimates in the wavelength range 


probed by the gamma-ray data. In order to account 
for the correlations between the gamma-ray based EBL 
points as well as for the broadness of the wavelength 
range including multiple points from galaxy counts, we 
build the following likelihood for the EBL excess, 
marginalizing over the EBL parameters and taking into 
account multiple galaxy counts estimates for a single 
gamma-ray-based EBL point: 




= J dA e 

^-^[A-Aoj-^VXdA-Ao] 


V % J 


X e 


(24) 


where Aq is the vector of the best-fit EBL parameters and 
where Vaq is their covariance matrix, provided for refer¬ 
ence in Appendix M and vlh are the EBL spe¬ 

cific intensity basea on gamma rays and galaxy counts, 
respectively. The integration is performed numerically 
through a Monte Carlo probe of the parameter space 
within ±3cr of the best-fit Aq. 

We measure an overall excess Ai/I^, = —0.8 ± 

0.4 nW m“^ sr“^, with a mild significance at the 2.0 cr 
level based on a likelihood ratio test. The slightly low 
gamma-ray estimate of the EBL is mostly due to the 
wavelength range above 25 /im, with values in the same 
units of 4.7 ± 2.2, 1.4 ± 1.4, -0.4 ± 0.6, and -0.8 ± 0.3, 
in the wavelength ranges 0.26 — 1.2 /im, 1.2 — 5.2 /im, 
5.2 — 23 /im, and 23 — 105 /im, respectively, as shown in 
Fig.S The coarse binning of the EBL spectrum based 
only on gamma-ray data is possibly responsible for the 
mildly negative excess, as no tension is observed in the 
eight-point spectrum when taking into account the local 
constraints. We note that the model developed for the 
excess of anisotropy at 1.1 /im and 1.6 / im detected with 


GIBE R by Zemcov et al. (2014) (see also Kashlinsky et al. 
2015 for a critical review), explained as originating from 


stars stripped from their parent galaxies during mergers, 
cannot be ruled out by our results given the uncertain¬ 
ties from both approaches. Little room is left for con¬ 
tributions from other unknown populations of sources, 
especially above 1.2 /im. This wavelength range is of par¬ 
ticular interest for the study of the sources of reioniza¬ 
tion and we exclude significant contributions fro m mini¬ 
quasars and Pop III stars as presented e.g. in jCoora^ 
& Yoshida (20041. We also exclude more exotic scenar- 
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Fig. 4.— Difference between the best-fit EBL spectrum derived 
from gamma-ray spectra only and the EBL estimates based on 
galaxy counts (blue squares). The anisotropy excess mea sured at 
1.1 Mm and 1.6/im and a model of the IHL (both from |Zemcov| 
|et al.||2014| l are shown with red-empty and dark-red fillea points, 
respectively. 


ios where Pop. Ill stars would experience a long “Dark 
St ar” phase powered b y WIMPs (see in particular Fig. 3 


Maurer et al.||20l2 ). 


4.2. Hubble constant 

The possibility to constrain the expansion rate of the 
universe using gamm a-ray observation of di stant sources 
was first proposed by Salamon et al. (1994). The idea is 
rather simple: measurements of the DDL optical depth, 
r, are proportional to vly/H^^ as shown in Eq. ^ while 
direct observations provide estimates of the EBL inten¬ 
sity vlv, the combination of which thus constrains Hq. 

Such an a pproach has been pursued e.g. by |Barrau| 
et al. (2008) who fixed the EBL intensity within esti¬ 
mates based on galaxy counts. Following a similar ap¬ 
proach, we assume here that the EBL intensity can be 
described by the lower limits in T able [l] This assumption 
follows the results from Sec. |4.1[ showing a rather good 
agreement between gamma-ray data and galaxy counts. 

We hereafter use solely gamma-ray data to estimate 
vI^/Hq. Calling H the true value of the Hubble con¬ 
stant and Hq = 70 km s“^ Mpc“^ the value used to de¬ 
rive the EBL parameters, Aq, we define the marginalized 
likelihood over the EBL parameters: 


C{H) = y 4^ 




X e 


(25) 


where XEBhi^) assesses the compatibility of the set of 
parameters A with l ocal constraints on the EBL, as de¬ 
scribed in Sec. 13.3.11 

The marginalized likelihood shown in Fig. yields 
an estimate of Hq = 88 ± 8stat ± 13sys km s“^ Mpc“^. 
The systematic uncertainty is propagated from the op¬ 
tical depth (7 — 8%, see Appendix 1^, with Hq oc 
1/t, and then added in quadrature ro the bias ex¬ 
pected from an excess of intensity with respect to galaxy 
counts . A n offset of —0.8nW m“^ sr“^, as determined 
in Sec|4.1[ yields an estimate of the Hubble constant 



Fig. 5.— Likelihood distribution of the Hubble constant. The 
estimate based on the comparison of gamma-ray data and local 
constraints is shown in solid blue. For reference, this estimate is 
compared to the best-fit values based on CMB observations and 
the local cosmic distance ladder. 


11 km s“^ Mpc“^ lower than the best-fit value. Account¬ 
ing for both statistical and systematics uncertainties on 
this measurement, no tension larger than 1.4 ct is ob¬ 
served with constraints based on the cosmic distance 
ladder (|Efstathiou||2014|) or CMB -based measurements 
(Planck Collaboration et al. 2014). 


Dominguez & Prada (|2U13I noticed that not only 
the distance term but also the density of photons used 
in the optical depth computation could depend on the 
Hubble constant for a given set of galaxy observations. 
Within such a formalism, the evolution parameter fevoi 
would be a function of Hq, with a typical span of 
0.5 < ffivni < 2.5 for 0.5 < Flg/lOOkm s“^ Mpc“^ < 0.9 
(Dommguez|2015 ). Marginalizing the likelihood over this 
range of evolution parameters yields mildly larger un¬ 
certainties {Hq = 88 ± 13stat ± 13sys km s“^ Mpc“^), 
but does not affect our conclusions. We note that 
the result of [Dominguez fc Prada (2013), Hq = 
71 -1-4.6 -1-7.2 Vm 0-1 Mpc“^, exploiting the cos¬ 


yi -I-4.D +I.-Z ^ 

— 5.6(stat) —13.8(sys) 


mic gamma-ray horizon (Blanch & Martinez 2005) and 
a fixed EBL spectral shape at z = 0 remains the most 
competitive gamma-ray estimate of the Hubble constant. 


4.3. Source redshifts 


A similar approach to that devised in Sec. 4.2 can 
be used to constrain the distance of unknown redshift 
sources. In the following, we use gamma-ray absorption 
as a guide for selecting possibly conflicting spectroscopic 
estimates. We assume a Hubble constant fixed to its 
nominal value Hq = 70km s“^ Mpc“^ and we describe 
the E BL b y its best-fit eight-point spectrum obtained 
in Sec. 4.1 (combined gamma-ray and local constraints). 
The best-ht parameters Aq and covariance matrix 
for the eight-parameter spectrum are provided in Ap¬ 
pendix [B| 

We marginalize over the EBL parameters within un¬ 
certainties and define the following likelihood for a given 
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Fig. 6. — Likelihood distribution of the redshift of six TeV 
blazars. Constraints are obtained after marginalization over the 
best-fit EBL parameters. 


set of spectra from a single source: 

L{z) = JdA Q-hyA-Ao^^vx^lA-Ao] 

^ g 2 points ~t XhE — VHE ) (^’^) ^26^ 

where the sum is over the spectra of a single source. Note 
that the spectra used in this sum must be different from 
the ones used in the estimation of the EBL parameters 
Aq, to avoid double counting the same datasets. The 
unknown-redshift sources are not included in the gamma- 
ray cosmology sample used to determine Aq, which jus¬ 
tifies our approach. All the intrinsic spectra studied in 
this section are well matched by intrinsic PWL models. 

Six sources with underconstrained distances are listed 
in Table [2] The likelihood distributions of the redshifts 
for the six sources are shown in Fig. 


lES 1215-1-303: for which we have gathered 
two spectra from MAGIC and VERITAS. Two 
spectroscopic estimates of the redshift of this HBL 
can be fou nd in the literat ure: z = 0.130 ( Bad e 
et al. 19981) and 2 = 0.237 (Lanzetta et al. 1993 T 


Akiyama et al. (2003) show a spectrum and list 
the redshift as 0.130, but it is not clear whether 
that redshift is supported by the spectrum or 
taken from the literature. As shown in Fig. the 
likelihood profile for this source has a maximum 
at z ~ 0.2. Because this estimate is compatible 
with zero, we provide an upper limit at the 99% 
confidence level of z < 0.35, and use in the 
following z = 0.237. 


PKS 0447-439: for which we have gathered one 
spectrum from H.E.S.S. A spectroscopic e s timat e 
of z = 0.205 was claimed by Perlman et al.| (1998), 
based though on a rather weak feature that was 
not confirmed b y furt her optical obser v ations 

(IMS) 


(e.g. fPita et al. 
recently estimate 
to be z = 0.343 


r2014|). 
a the 


Muriel 
redshift 




et al. 
this 


source 

based on the observation of 
neighbouring galaxies possibly belonging to the 


same cluster. The H.E.S.S. spectrum has been 
used by several teams to constrain the distance 
of this source, with claimed me asurements of 
z = 0 . 16±0.05 arid z = 0.20 ±0.05 (|Prandini et al.j 
TTVToi ir7i — Horn Al respectivel;^! -n 


2012 iZhou et al. 2014 


Such small 

uncertainties, related to limiting assumptions on 
the EBL model or on the properties of the intrinsic 
spectra of TeV blazars, are not confirmed by the 


H.E.S.S. Collaboration 
limit 


(2013b), 


which provides 
level of z < 0.59. 


an upper limit at the 
Our profile for PKS 0447-439 peaks at z ~ 0.2 
and we obtain z < 0.45 at the 99% confidence 
level. In the following, we use z = 0.343 for this 
source, noting that z = 0.205 yields compatible 
results given the broadness of the likelihood profile. 


S5 0716-1-714: for which we have gather e d one 


spectrum from MAGIC. Danforth et al. (2013) 


detected a system on the line of sight and have 
also set an upper limit based on the non-detection 
of lines farther away, constraining the redshift 
of this IBL within 0.232 < z < 0.322. This 
range is consistent with the spectroscopy of three 
galaxies possibly hosted by the same cluster as th e 
source, around z ~ 0.26 (Bychkova et al. 2006). 
In our study, no conclusion can be drawn tor 
S5 0716-1-714, which shows a very broad maximum 
around z ~ 0.5. In the following, we use the 
cluster-based z = 0.26 as the redshift estimate of 
this source. 


3C 66A: for whi ch we have gath e red one spectrum 


from MAGIC, 
an absorber 


on 


Furniss 
the 


line 


et al. 

3F 


( |2013a ) 
sight and 


detected 
have con¬ 


strained the maximum distance of this IBL based 
on the non-detect ion of farther lines, finding 
0.335 < z < 0.41. |Yang & Wang| 12010) studied 
the gamma-ray emission of the object and set 
an upper limit at z < 0.58. For 3C 66A, the 
likelihood in Fig. peaks at a rather low redshift 
(z ~ 0) and indicates a distance fully compatible 
with zero within uncertainties. No upper limit on 
the redshift of 3C 66A can be obtained given the 
broadness of the distribution. In the following, 
we use the spectroscopic lower-limit of z = 0.335, 
which is only in mild tension with the likelihood 
profile (at the 1 cr level) 


• PG 1553-1-113: for which we have gathered two 
spectra from H . E.S.S . and three from MAGIC. 


Danforth et al. (2010) constrained the redshift of 
this HBL through spectroscopic observations to 
0.433 < z < 0.58, though the lower limit is based 
on a single line. The authors obtained z > 0.395 
based on more numerous absorbers. Ba s ed on 
gamma-ray observations, |Yang & Wang| (|2010| 
determined z < 0 .78. Fixing the EBL absor ption 


to the model of Franceschini et al. _([2008 ), the 


H.E.S.S. collaboration 
measured z = 0.49 ± 


Abramowski et al. 
.04, or z 


< 


2015) 
0.56 at the 


95 % confidence level, compared to z < 0.60 and 
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z < 0.62 from the MAGIC and VERITAS collabo- 



yields 

conservative gamma-ray estimate of the redshift of 
this source, with most likely value z = 0.4lto;?i, 
preferred at the 3.4 cr level to a null redshift. 


For compari son with the work of Abramowski 


et al. (2015), we obtain a 95% upper-limit on 


the redshiit of the source z < 0.53 (z < 0.58 at 
99%). In the following, we set the redshift of PG 
1553-fll3 to z = 0.433. 


PKS 1424-1-240: for which we have gathered three 


Furniss et al. 

(2013b 

) determined a spectroscopic 

lower limit on the redshiit ol this object, z > 0.604. 

Prandini et al. ( 

2011 

1 claim a measurement of 


from the s ame limiting assumpt ions as for PKS 
0447-439. Yang & Wang (2010) provide a more 
robust upper-limit at z < 1.19. As for 3C 66A, the 
likelihood proHle for PKS 1424-1-240 is fully com¬ 
patible with zero redshift, but an upper limit can 
be set at z < 0.64 at the 99% confidence level. In 
the following, we use the spectroscopic lower-limit 
of z = 0.604, which is in slight tension (2.4 cr) with 
the likelihood profile. 

The four constraints on redshifts derived in this section 
are the most stringent gamma-ray upper limits to date 
for these sources. 

4.4. Axion-like particles 

With a best-fit EBL spectrum at hand and redshift 
estimates for all sources, we can search for any siginih- 
cant residual, which might be indicative of new physics. 
Though we restricted the study to a gamma-ray cosmol¬ 
ogy sample of spectra for the determination of the EBL 
and of the Hubble constant, we study here all the spec¬ 
tra listed in Table with spectral models provided in 
column 5. 

We investigate in this section deviations from classical 
absorption resulting from coupling of gamma rays with 
axion-like particles. Gamma rays from blazars could con¬ 
vert into such hypothetical particles (cousins of the QCD 
axion with free mass and coupling) and could then con¬ 
vert back into gamma rays within the Galactic magnetic 
field. Several locations have been studied for the initial 
conversion, be it within the magnetic field of the emis¬ 
sion region, in the vicinity of the source (host, jet galaxy, 


3 


cluster) or in the intergalactic medium ( see e.g. Hooper 
&: Serpicol 2007[ |de Angelis et al.ir2007| |Sancfiez- Gonde 
et al. 112009 ITavecchio et al.||2012D. An observabi 


Dservable effect 
would be a gamma-ray spectral hardening at the highest 


eral analyses (Dominguez et al. 2011a Horns & Meyer 

2012 

iVleyer et al. 

2013 Rubtsov & 'lroitsky||2014). We 


note nonetheless that a large traction of the parameter 
sp ace corresponding to thes e hints has been excluded 
by Abramowski et al. (2013) using the lack of point-to- 
point fluctuations in the spectrum of the bright blazar 
PKS 2155-304. 



10 ' 1 10 


Optical depth x 

Fig. 7.— Residuals to the best-fit models for the 106 spectra (737 
points) studied in this paper, as a function of optical depth. 


The spectral points and best-ht models are shown in 
Appendix [C| We show in Eig. [^the residuals to the best- 
fit models as a function of optical depth for these 737 
points. All of the spectra are well represented by our 
models. The residuals are scattered around an average 
of 0.05 ± 0.03 with an rms of 0.73 ± 0.02 and are well 
represented by a normal distribution, with a p value of 
22 % following an Anderson-Darling test. The rms is 
signihcantly smaller t han 1 . which is consistent with the 


small found in Sec. 4.1 


In the following, we search for a harde ning of the TeV 
spectr a at the highest optical depths as in [Horns fc Meyer] 
(20121, while leaving model-dependent constraints on the 


axion-lik e particles’ coupl i ng to gamma rays to future 
studies. Horns & Meyer (2012) model a sample of 50 
sp ectra and assume an EB L intensity fixed to the model 
of Kneiske & Dole (2010). Noting an rms of the flux 
residuals smaller than one, and deviations from a normal 
distribution b ased on an Anderson-Darling test. Horns & 
Meyer (2012) chose to ignore uncertainties on the data 
and to study the following quantity: 


R = 


4^i 0 model ,2 
T 0model,2 


(27) 


where (j)i is the observed flux for the point i and (/)modei,i 
is the expected flux resulting from the model. 

Gomparing the distribution of i? in a reference sample 
composed of points at 1 < r < 2 and in a search sample 
at r > 2, the authors find a 4.2 a discrepancy, with an av¬ 
erage value of i? ^ 0.25 at large optical depth (estimated 
from the cumulative distribution function shown in Eig. 3 
of their publication). Dehning the flux-enhancement fac¬ 
tor as FE = (j)/(jimodeh the flux of TeV blazars would then 
exceed the expectations by FE = (l-|-i?)/(l — i?) ~ 1.7 
above an optical depth of r > 2. They explain this 
^ 70% increase in flux as an hint for mixing of gamma 
rays with axion-like particles. Performing the same test 
with our best-fit EBL spectrum and our larger dataset, 
we find a slightly larger discrepancy of 4.5 tr, for an av¬ 
erage R ~ 0.1. However, we argue in the following that 
this test is flawed, because it neglects the uncertainties 
on the measured flux. 

We measure the average flux enhancement EE in var- 
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Fig. 8 .— Flux enhancement, defined by the ratio of observed 
and expected fluxes, as a function optical depth. The sh aded gray 
reffion is the flux enhancement implied by the results oflHorns &| 


ious optical-depth bins, using a sample of gamma-ray 
spectra twice as large as any other studied. The residuals 
show a normal distribution, whose parameters are largely 
independent of redshift, energy, and optical depth. We 
thus include the uncertainties in the computation of the 
flux enhancement, weighting the relative contributions of 
the spectral points, as: 


FE = 




^model 




^model, 


/c 


(28) 


which is the usual y^-based weighted average, propagat¬ 
ing the uncertainty on the observed flux Note that 
scaling the uncertainties on the observed flux up or down 
affects the errors on the flux enhancement estimate, but 
not the mean. 

The flux enhancement is shown as a function of op¬ 
tical depth in Fig. The 106 spectra constrain the 
flux enhancement at optical depths larger than r > 2 to 
FE = 0.98 ± 0.04, for an average optical depth t = 3.0. 
This value, consistent with 1, shows no deviation from 
expectations. A flux enhancement of more than ^ 40 % 
is excluded beyond the 5 a level, taking only statistical 
uncertainties into account. Even assuming a systematic 
bias of 10% at large optical depth (see Appendix [^, the 
flux enhanceme nt corresponding to the results o: 


orns 


& Meyer| (|2012|) (gray shaded region in Fig. remains 
excluded at the 5 a level. 


Applying the statistical test of Horns & Meyer (20121 
to our dataset with our EBL spectrum therefore shows 
a more significant effect (4.5 ct compared to 4.2), albeit 
with smaller amplitude than they observe (i? ~ 0.1 com¬ 
pared to ~ 0.25). However, we find little indication that 
the uncertainties on the individual flux measurements are 
unreliable, and when they are taken into account, we do 
not find a flux enhancement at large optical depths. The 
highest flux enhancements at optical depth above 2 in the 
sample studied here are obtained for the last points of the 
spectra of lES 1959-b650 (Whipple, 2002), at r = 5.2, 


3C 279 (MAGIC, 2006), at r = 4.2, and lES 0229-b200 
(H.E.S.S., 2005-2006), at r = 6.7. They show flux en¬ 
hancement of 13, 14 and even 23, respectively, but with 
uncertainties o n the order of 100 % tha t strongly bias the 


test devised by Horns & Meyer (20121. The most likely 
source of anomaly seen by these authors seems then to be 
the choice of statistical test rather than a specific dataset 
(their set of spectra is to a large extent inc luded in Ta¬ 
ble |2l) . We note that |Meyer et al.| (|2013|) studied the 
residuals normalized to uncertainties, as shown in Fig.l^ 
and found an average flux enhancement significant at the 


4.4 CT level based on a t test and using the EBL model of 

Kneiske & Dole 

(2010 

). Performing the same test with 

our sample vield 

s a simiticance of 1.1a for the model of 

Kneiske & Dole 

(^2010 

) and 2.6 CT for the EBL spectrum 
the slightly different response of 

derived in Sec. 

4.1| ' 


due to the weighting of the points in the averaging pro¬ 
cess (oc 1/cr^^i for the t test and oc l/c^ j for a y^-based 


average in Eq. 28). 

We conclude, based on the largest gamma-ray sample 
studied so far, that current VHE gamma-ray observa¬ 
tions do not show a detectable flux enhancement at large 
optical depths and find little motivation for a lower limit 
on the cou pling of axion-like p articles with photons as 


reported in Meyer et al. (2013). 


4.5. Lorentz invariance violation 

We investigate in this section the compatibility of the 
86 spectra in the gamma-ray cosmology sample with a 


fune ( 

1999): Aloisio et al. 

20001); Protheroe & Meyer 

(2000 

; Ellis et al.l (|200l|): 

Amelino-Gamelia & Piran 

(2001) 

; Stecker & Glashow 

(iJOOl). These authors con- 


the s pectrum of Markarian 501 observed by HEGRA in 
1997 ( Aharonian et al.|1999 ), but the limited constraints 


on the EBL and on blazars’ spectra, which were only 
a handful in 2000, prevented a quantitative constraint. 
This spectrum is not included in our study because it 


was updated in Aharonian et al. (2001). 

We model the effect of a Lorentz invariance violation 


(LIV) adopting the formalism ofjJacob & Piran (2008). 
The starting point consists in a leading-order niodihca- 
tion of the special relativistic relation 
where E and p are the energy and momentum of a par¬ 
ticle of mass m. The effect should become significant 
around a quantum-gravity energy scale Eqq, a correc¬ 
tion of order n = 1, 2 yielding: 


E^ =p^c^+m^c'^±E^ ( (29) 
\EqgJ 

so that the norm of the momentum four-vector, E‘^—p^(?, 
is not a Lorentz invariant any more. 

Equation alters the kinematics of electron-positron 
pair creation, as in the collision of a TeV gamma ray in¬ 
teracting with an EBL photon. One needs to modify the 
special-relativistic threshold ethr > mgC^/E, with rUe the 
mass of the electron and where e and E are the energy 
of the two photons. The original spectrum of Markar¬ 
ian 501 published by HEGRA hinted at a gamma-ray 
absorption lower than predicted by classical interaction 
with the EBL, which corresponds to the subluminal case 
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correction in Eq. 29) that we study in the following. 
A superluminal correction (“+” in Eq. 29) would yield 


either a lower thresho ld and pair creation on the CMB 
(Jacob & Piran 2 008) or photon de cay if the dynamics 
are favorable (e.g. |bhao & Ma||2010|), in which cases the 
absorption would be larger than classical. The investi¬ 
gation of the superluminal scenario, possibly requiring 
assumptions on the dynamics of the process, is left for 
future studies. 

In the subluminal case, the energy-momentum conser¬ 
vation yields the modified pair-creation threshold: 


2 4 

mtc 


1 - 2 - 


^thr — 


E 




(30) 


if Eq.|^is applied to the photons and leptons. The term 
1 — 2“” should be replaced by 1 if only the photons are 
affected, which is equivalent to considering a quantum- 
gravity energy scale twice as large for n = 1. One can 
rewrite Eq. [30|as: 


^thr — 


2 4 

mtc 


E 


E 


E, 


7,LIV 


n+2 


(31) 


where 


n\l 1/("+2) 


/ E \ 

= 29.4TeV x , for n = 1 (32) 

\ -^Planck / 


with Epianck = jG = 1.22 X 10^® eV. A leading 
quadratic correction, n = 2, yields E^^isi = 120 PeV x 
a//T-P ianck, out of reach of current experiments for 
Eqg ~ Epianck- We focus in the following on the linear 
case, n = 1, but we derive the modified EBL optical 
depth in the general case as: 


t{Eo,zo) 


3 CTpC 






de^(e,z) 



(33) 


where we follow the notations of Sec. [21 
Assuming that the energy and redshift dependences of 
the EBL photon held can be decoupled and using the 
same changes of variable as in Sec. the LIV-affected 
EBL optical depth can be expressed as: 


_ , 37raT Eo , 

t(Eo,zo) =—— X X / deo exp(-3eo) 

-^0 Tfl^C J QQ 


X uly ( eo — In 


En 


dl evol(z) 
dz- 


J Jo dz(\ + zY 


1-^ 


(1 -)- z)Eo 


E. 


where 


2 exp(-eo) 

Pmax ^^2 


7.LIV 


1-h 


n+2 


PiPn 


(1 -)- z)Eo 


E. 


7,LIV 


n+2 


(34) 


(35) 


We compare the LIV-affected optical depth to the clas¬ 
sical one in Fig.|^ using the eight-point specihc intensity 



Fig. 9.— EBL optical depth at three redshifts, in the classical 
case (solid lines) and in the case of a LIV modification at the Planck 
scale (dashed lines). 


of the EBL obtained in Sec. 14.11 The main effect of LIV 
on gamma-ray absorption is a reduction of the optical 
depth above lOTeV, largely indep end ent of redshift. 

The approach developed in Sec. |3.3| consists in hnding 
the best EBL spectrum jointly accounting for the ab¬ 
sorption signature seen in gamma rays and for local con¬ 
straints on the EBL. We add here an extra free parameter 
that is the energy scale Eqq at which LIV modifications 
of the pair-creation threshold take place. We vary Eqq 
and compute the best-fit accounting for both gamma- 
ray spectra and local constraints on the EBL. We de¬ 
fine the test statistic TS = (Eqg, Aqq ) — (oo, A^o ), 
where the latter x^ is measured at the best-fit EBL level 
in the classical case, and where A denotes the EBL pa¬ 
rameters, left free to vary for each quantum-gravity scale. 

Figure shows the likelihood profile, C = 
exp(—TS'/2), as a function of the inverse of Eqq normal¬ 
ized to the Planck energy. Interestingly, a slight excess 
can be seen around the the Planck scale, though with a 
significance of only 2.4 cr. Our study is performed using 
a sample of 86 spectra that includes the very observa¬ 
tion of Markarian 501 by HEGRA that triggered sub¬ 
stantial theoretical developments on modifications of the 
pair-creation threshold. The original spectrum, which 


Aharonian et al. 

(1999 

). The data were subsequently re- 

analyzed by Aharonian et ai. 

(2001) with an improved 

energy resolution, yielding a compatible spectrum cover- 


ing the energy range 3.3— 21 TeV, which is the spectrum 
included in our gamma-ray cosmology sample. Using in¬ 
stead the original spectrum enhances the Planck-scale 
excess to the 4 ct level. The smaller energy coverage 
of the higher-resolution spectrum results in a less con¬ 
strained intrinsic emission and the observed spectrum 
has a somewhat sharper rollover at the highest energies, 
explaining the decrease in significance from 4 ct down to 
2.4 tr. Given the low significance of the excess, we pro¬ 
vide lower limits on the quantum-gravity energy scale of 


Eqq > 0.78 X Epianck (95%) and Eqq > 0.65 X Epianck 
(99%). A 10% systematic uncertainty on the energy 
scale of gamma-ray instruments similarly impacts the 
quantum-gravity energy estimate, and we quote a robust 
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Fig. 10.— Likelihood profile of the quantum-gravity energy scale, 
leaving the eight-point EBL spectrum free. 


lower limit accounting for the systematic uncertainties of 
EqG > 0.6 X iipianck. 

Other constraints on LIV have been derived from ob¬ 
servations of high-energy sources. The synchrotron emis¬ 
sion by electrons in the Crab n ebula has been used, 
e.g. by Jacobson et al. (2003), to constrain a lin¬ 
ear subluminal term (the same as constrained here) 
to a quantum gravity scale seven orders of magnitude 
above the Planck scale. Some authors nonetheless ar¬ 
gue that this test is not only kinematics but also de¬ 
pends on the dynamics of process, which are far from 
bein g understood within qu antum-gravity phenomenol¬ 
ogy ( Amelino-Camelia|2013 ). A suppressed cross section 
would then mitigate the constraints from the Crab neb¬ 
ula. A more pristine test comes from observations of 
gamma-ray bursts (GRBs) by Fermi, assuming that LIV 
affects the propagation time of photons in an energy- 
dependent way. The most recent constraints have been 
derived by Vasileiou et al. (2013) who obtain, after tak¬ 
ing into account possible internal delays in the source, 
a lower limit Aqg > 2 x Apianck from a single GRB 
(090510), and Age > 0.1 x Apianck from three others. 

One could conclude from the GRB 090510 limit than 
any modification of the pair-creation threshold is ruled 
out up to 2 X Apianck, but this would be overlooking the 
fact that time-delay and absorption observations con¬ 
strain two different processes. The underlying theory 
might indeed preserve the speed of light, not impact¬ 
ing time delay observations, but at the same time affect 
the dispersion relation of particles, leaving an imprint 
in blazar s’ spectra (i.e. v = dE/ dp would no longer be 
valid, see Amelino- Camelia||20 13). High-statistics obser¬ 
vations above lU'l’eV will be required to further study 
quantum-gravity effects, and the formalism we have de¬ 
veloped in this section will prove useful to further con¬ 
strain LIV with gamma-ray observations of blazars. 


5. CONGLUSIONS 

A wealth of data has been published by ground-based 
gamma-ray observatories over the past two decades. We 
have compiled the most extensive set of gamma-ray spec¬ 
tra from VHE blazars to date, with 106 spectra from 38 
objects, corresponding to a total of about 300,000 gamma 


rays. This is twice the size of any sample studied before. 

Our first result is purely analytical. We have discov¬ 
ered that the triple integral relating the gamma-ray op¬ 
tical depth to the EBL intensity can be reduced to a 
double integral without any loss of generality. We have 
further shown that, assuming a decoupling of the evo¬ 
lution of the EBL and of its rest-frame spectrum, the 
gamma-ray optical depth is the convolution of the EBL 
intensity with the EBL kernel. 

This analytical work significantly reduces the complex¬ 
ity of the spectral reconstruction of the EBL based on 
gamma-ray observations. The decoupling approxima¬ 
tion introduces rather small (< 5 %) systematic errors for 
sources in the local universe. Using a joint spectral analy¬ 
sis of a subsample of 86 spectra, we deconvolve the intrin¬ 
sic emission of the sources from the imprint of the EBL 
spectrum. The reconstructed EBL intensity is preferred 
at the 11 (T level to the absence of gamma-ray absorption, 
and we reconstruct an eight-point spectrum covering the 
wavelength range 0.26 — 105 /rm, from mid-UV to far IR. 
The spectrum of the EBL based on gamma-ray obser¬ 
vations is in good agreement with estimates based on 
galaxy counts, with uncertainties that lea ve some room 


for co ntributions from e.g. intra-halo light (Zemcov et al. 


2014), while constraining the emission of reionization 
sources such as Pop. Ill stars or miniquasars. The 
brightnesses of the COB and GIB are measured to be 
36 ± 11 nW m“^ sr“^ and 25.9 ± 3.4nW m“^ sr“^, re¬ 
spectively. Once integrated between 0.1 — 1000 ^m, the 
EBL is 6.5 % of the CMB, with an overall uncertainty on 
this number of 20 %. 

Our third result is a gamma-ray measurement of the 
Hubble constant, iLo = 88±8stat±13sys km s“^ Mpc“^, 
that is both model independent and based on a signih- 
cant number of gamma-ray spectra. Such constraints on 
the expansion rate of the universe are independent from 
measurements based on the CMB and the cosmic ladder 
of distances. 

We measure no significant flux enhancement at large 
optical depths and we rule out at the 5 a level the “pair- 
produ ction anomaly” as obtained by [Horns fc Meyer] 
(2012). This suggests that the level of mixing of axion- 
like particles with TeV photons, if any, is below previous 
estimates based on this effect. We would also like to cor¬ 
rect two misconceptions sometimes circulating in the lit¬ 
erature. (i) The best-fit EBL spectrum based on gamma- 
ray observations is not significantly lower than the min¬ 
imum EBL level, (ii) The intrinsic gamma-ray spectra 
reconstructed after deabsorption from the EBL effect are 
not too hard with respect to expectations. All the intrin¬ 
sic spectra in our fits are softer than measured at lower 
energy when contemporaneous data are available. This 
diminishes the motivation for scenarios where axion-like 
particles impact gamma-ray absorption, or where secon¬ 
daries from UHECR contribute to the gamma-ray signal 
observed on Earth. 

Einally, we constrain the impact of a Lorentz invari¬ 
ance violation on gamma-ray absorption. A modified 
dispersion relation, with a correction scaling as the ratio 
of the gamma-ray energy and the Planck energy, alters 
the threshold of pair creation and results in milder ab¬ 
sorption of > 10 TeV gamma rays from blazars. Our 
formalism takes into account both the spectrum and the 
evolution of the EBL. It also enables a quantitative anal- 
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ysis of blazars’ spectra in the presence of LIV. A weak 
2.4 cr excess prevents us from ruling out a modification at 
the Planck energy, but we rule out for the hrst time any 
effect below 0.6 x Apianck at the 99 % confidence level. 


Th e successes of gamma-ray cosmology (see Biteau 
2013 for a review) are far from being closed by this work. 


Major progress in UV and IR observations of distant 
galaxies has been achieved in recent years, in particular 
with Galex, Spitzer, and Herschel. Further achievements 
are exp ected from the Jam es Webb Space Telescope 
(JWST,|Gardner et al.|2006|) during the next decade be¬ 
tween 0.6 and 27/rm, with improved constraints from 
galaxy counts in the optical and mid-IR. The combina¬ 
tion of next-generation local constraints with the tremen¬ 
dous gamma- ray sensitivity of the Chenrekov Telescope 
Array (CTA, |Acharya et al.||2013fc above 30 GeV could 
signihcantly rehne the binning of the EBL spectrum. 
The theoretical limit is dictated by the energy resolution 
of gamma-ray telescopes, about 10 % in gamma-ray en¬ 
ergy or equivalently 10 % in EBL wavelength. Such fine 
spectroscopy could clearly show spectral signatures from 
polycyclic aromatic hydrocarbons, from intra-halo light, 
or from the sources of the reionization of the universe. 
The combined constraints of GTA and JWST will also 
be crucial for the measurement of the Hubble constant 
based on gamma-ray absorption. 

Further developments of analyses at the event level 
(binned or unbin ned), such as 3ML or GammaLib 
( Knodlseder|2013 ), combined w ith full likelihood s pectral 


techniques such as developed in 
increase the statistical power o: 


Piron et al. (20011, could 
the approach developed 


in this paper. Such methods are in principle not affected 
by overestimated uncertainties from published spectra, 
where correlations are neglected, and would open the 
possibility of a broad-band spectral fit, e.g. combin¬ 
ing the data from fermz-LAT, G herenkov telescopes, and 
HAWG (Abeysekara et al. 2013). The subheld of gamma- 
ray cosniology focusing oii the rate of the pairs produced 
by gamma-ray absorption, and how they cou ld be im¬ 
pacted by the intergalactic mag netic held (see Durrer & 
Neronov 2013 Chen et al. 2014), will greatly beneht from 
these ongoing developments. Finally, we are working on 
the extension of the method that we have presented here 
to higher redshifts (z > 1), where a simple parametriza- 
tion of the EBL evolution starts to fail. If successful, sig- 
nihcant improvements in the spectrum of the EBL below 
1 /im can be expected from the analysis of Eer77iz-LAT 
blazars. 
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Fig. 11.— A ttenuation curves for thr ee differen t redshifts, from 2 = 0.03 up to 2 : = 0.6. Gray curves are directly extracted from the 
publications of|Franceschini et al.|(|2008|| (left) and|Gilmore et al.|(|2012|| (right). Colored curves are based on the same z = 0 EBL density 
as in the publications but assume a template EBL evolution witn fevoi = 


APPENDIX 

A. SYSTEMATIC UNCERTAINTIES 


We discuss in Sec. |A.1| and S ec. |A.2| the approximations to the spectrum and evolution of the EBL used in this 


paper. We estimate in Sec. A.3 the systematic uncertainty arising from these approximations as well as those coming 


from the modeling of the intrinsic spectra and the possible biases in the energy scale. 

A.l. EBL evolution 


We compare, in Eig. 11 gamma-ray optical depths published by Franceschini et al. (2008) and Gilmore et al. (2012) 


with optical depths derived from Eq. assuming the z = 0 specific intensities as given m each paper with a ternplate 
evolution with /evoi = 1-7- The publisned attenuation curves and those derived from Eq.J^are in very good agreement, 
which supports th e an alytical approach in Sec.The differences between published ancftemplate optical depths. At, 


are shown in Fig. 12 The value of /evoi = 1-7 is adopted (intermediate panels) and can be c ompared to sof t er an d 
harder evolution imtne bottom and top panels, respectively. We note that the evolution used by Raue & Mazin| (2008), 
/evoi = 1-2, results in significant deviations at large redshifts with respect to the models. Similarly, the published optical 
depths are underestimated by the template approach at large redshifts for a soft evolution with /evoi = 2.2. 

The template evolution with /evoi = 1-7 results in an optical dep th difference on the order of 0. 1 on average, which is 
comparable with the difference in evolution between the models of Franceschini et al. (2008) and Gilmore et al. (2012) 
themselves. As far as EBL evolution is concerned, assuming an energy-redshitt decoupling in the local universe has 
then a similar impact on the absorption to using one or another state-of-the-art model. 

Below an optical depth of 3, the deviation in the absorption factor remains smaller than 15%, which is below the 
typical systematic errors on gamma-ray fluxes measured by current-generation ground-based instruments. Another 
reference point for the difference in optical depth can be provided no ting that a 5% deviation in Hn, r oughly the 


difference between Hq — 70 km s ^ Mpc ^ and the best-fit value from Planck Collaboration et al. (2014), results in 


a 5% deviation in optical depth, which corresponds to At = 0.15 (15% error on the absorption) for r = 3. Thus the 
template approach that we use in this publication introduces errors in the EBL no larger than those resulting from 
the uncertainties on Hq or from the differences between state-of-the-art models. 

For reference, the integrations in Eq.j^and Eq. |10[ which enable the computation of the optical depth, are performed 
using the trapezoidal rule, with uniforms steps in redshift of Sz = 10“^ and in logarithmic reduced photon energy, 
6eo = 5 X 10“^, from the threshold of the EBL kernel, eg = —21n(l -|- zq), up to eg = 10^. We checked that above 
eg > 10^ the tail of the kernel in Fig. has a negligible contribution to the convolution product in Eq. The 
overall uncertainty on the optical depth T(i?o,zg) due to the numerical integration is on the order of 0.01-0.03, mildly 
depending on the redshift of the gamma-ray source. 


A.2. Gaussian-sum approximation 

Besides the evolution of the EBL, our second sourc e of systematic error com es from the approximation of the true 
spectrum of the EBL by a sum of Gaussians, as in Eq. [T^ We show in Fig.[^the approximation of a “smooth” model 
vly{\), namely a sum of two log parabolas, by a Gaussian sum for a binning A; = 0.75, between 0.26 — 105 /im, as 
in Sec. and The weights of the Gaussian-sum model are obtained by numerically inverting Eq. |20[ where we set 
vll = 

The EBL kernel in Eq. [^smooths the EBL density over a wide range of EBL wavelengths. The small deviations in 
intensity arising from the uaussian-sum approximation then result in even milder optical-depth deviations, typically 
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Fig . 12.—- Left: Difference between optical depths obtained with a template evolution and with the model published by I Franceschini 
et al. l|2008l l. Right: Difference between optical depths obtained with a template evolution and with the model published by iGifmore 
et al. From top to bottom, the evolution parameter is /evoi = 1-2, /evoi = 1-7, and /evoi = 2.2. The intermediate evolution witn 

= 1.7 IS adopted. 


St < 0.1, 0.01, and 10 ^ for A; = 1.0, 0.5, and 0.1 respectively. 


A.3. Quantifying the systematic errors 

To estimate the systematic uncertainties arising from the modeling of the E BL, i.e. from the template evol ution and 
the Gaussian-sum approximation, we compare the optical depths derived by Franceschini et al. (|2008[) and Gilmore 


et al. (2012) to optical depths obtained with template evolution and approximating the Franceschini et al. (2008|) or 
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Fig. 13.— Smooth EBL intensity (sum of two log parabolas, solid cyan curve) approximated between 0.27 — 105 fira by a sum of Gaussians 
of width A; = 0.75 (dotted light blue). The resulting intensity is shown as a dark-blue dashed line. 


Gilmore et al. (2012) SEDs at z = 0 by Gaussian sums. 

We weight the contributions of the different optical depths based on the uncertainties on the spectral points included 
in the analysis. The measured optical depth depends on the measured flux, cj), as Ini/i = Int/iint —t, so that the maximum 
uncertainty on the optical depth, obtained by fixing the intrinsic model, scales as cr^/^. Then, one can estimate the 
EBL normalization factor a that accounts for the change from the model, of optical depth Tmodeh to the template 
approach, of optical depth rtempiate, by minimizing: 






(iteniplate,! H X Imodefi) 


which gives 


^model,i^template,i X {4‘il 


E^-LdeU X 


(Al) 


(A2) 


Since the optical depth is proportional to the EBL intensity and inversely proportional to the Hubble constant, the 
systematic relative bias of the optical depth, a — 1, contributes equally the systematic uncertainty of the EBL intensity 
and Hubble constant. Using the ga mma-ray cosmology sampl e, and assuming a binning of A/ = 0.75, we fin d an 
average bias of 0.5% for the model of Franceschini et al. (2008) and 1.8% for the model of Gilmore et al. (2012). 


TABLE 5 

Systematic uncertainties on the EBL specific intensity estimated 
WITH TWO different MODELS. 


Franceschini et al. 

2008 


Gilmore et al. 

2012 


EBL modeling 

0.5% 

1 .8% 

Intrinsic model 

2.3% 

5.2% 

Energy scale 

6 .2% 

6 .0% 

Total 

6 .6% 

8 .1% 


This source of systematic uncertainty is compared with the two other principal sources that have been identified: the 
gamma-ray energy scale and t he choice of intrinsic model for each spectrum . Th e EBL scaling factor of the models of 


Franceschini et al. (2008) and Gilmore et al. (2012) are estimated as in Sec. 4.1 additionally shifting the energy scale 
ot the Cherenkov experiment within ±10% on one hand, and using exclusively log-parabolas or exponential cut-off 
power laws for the intrinsic spectra on the othe r hand. A sy s temat ic energy bias of ±10% is a conservative estimate 
with respect to the value half as large found in Meyer et al. (2010) by studying the marginal mismatch between HE 
and VHE observations of the Crab nebula, but is comparable to the systematic error typically quoted by the observers. 
The energy scale impacts the EBL intensity at the ^ 6 % level. Changing the gamma-ray spectral models impacts the 
EBL normalization at the 2 — 5 % level. The total systematic uncertainty on the EBL flux level and on the Hubble 
constant are estimated in Table by summing in quadrature the various sources of systematic errors to be about 
7-8%. 

We also estimate the systematic bias on the flux enhancement discussed in Sec. 4.4 Using the full sample, we 


compute the average flux bias, FB, as the flux enhancement resulting from using the template approach instead of 
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Fig. 14.— Ratio of gamma-ray attenuation from the model of |Franceschini et al.| l|2008|) (light blue circles) andjGilmore et al.|(|2012j 
(dark blue squares) to the value obtained with our templates, using the lull sample ot gamma-ray spectra. 


the model: 


FB = exp[(a - 1) X Tmodel] 


(A3) 


where a is computed as in Eg. |A2[ The flux bias is shown as a function of optical depth in Fig. [T^ for the full sample. 
Using the gamma-ray cosmology sample yields similar average values for each of the two models, though with a rather 
constant behavior across the whole optical depth range. An underestim ation of the absorption of about f 0 % at optical 
depth r > 2 is visible in Fig. 14 for the model of Gilmore et al. (20121, although the large uncertainties do not allow 
a firm conclusion. For the sake of argument, we consider 10% as the systematic uncertainty on the flux enhancement 
at large optical depths. Improved statistics will result in a more accurate estimate of this systematic uncertainty in 
the future. 


B. BEST-FIT PARAMETERS AND COVARIANCE MATRICES 

In the following, we describe the best-fit EBL spectrum derived using gamma-ray data only (4 parameters) and 
gamma-ray data together with local constraint s (8 parameters). The EBL intensities and associated uncertainties in 
each wavelength bin can be computed from Eg. |20| The best-fit parameters of the gamma-ray data only are shown in 
Table with covariance matrix shown in Table 7 The best-fit parameters and covariance matrix for the gamma-ray 
data and local constraints are shown in Tables 8 and T he optical depths derived with the eight-point spectrum 
between 50GeV and 20TeV are shown in Tables|l0| and 1 1 1] We limit the results to optical depths less than 5, where 
we have good control of the systematic errors, as oTscussot in Appendix |A] 


TABLE 6 

Parameters of the best-fit EBL spectrum 
USING GAMMA-RAY DATA ONLY. 


A 

di 


nW sr“^ 

0.55 

18.1 

2.47 

8.05 

11.1 

0.91 

49.7 

4.51 
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TABLE 7 

Covariance matrix of the best-eit EBL parameters using gamma-ray 

DATA ONLY. 


1 2 

3 

4 

18.66 4.58 

1.62 

0.38 

1.83 

0.52 

0.17 


0.35 

0.04 



0.21 


TABLE 8 

Parameters of the best-eit EBL spegtrum using gamma-ray data and 
LOCAL constraints. 


A 

ai 

nW sr “^ 

0.38 

7.12 

0.80 

12.6 

1.70 

10.3 

3.60 

6.79 

7.62 

- 0.61 

16.1 

4.56 

34.1 

1.30 

72.3 

11.9 


TABLE 9 

Covariance matrix of the best-eit EBL parameters using gamma-ray 
DATA and local CONSTRAINTS. 


1 2 

3 

4 

5 

6 

7 

8 

21.710 - 3.450 

0.659 

- 0.001 

- 0.054 

- 0.007 

- 0.006 

- 0.020 

6.153 

- 0.435 

1.244 

0.076 

0.076 

- 0.001 

0.019 


0.802 

- 0.108 

0.232 

- 0.004 

0.012 

0.031 



0.827 

- 0.069 

0.066 

- 0.008 

0.004 




0.434 

- 0.035 

0.032 

0.052 





0.109 

- 0.030 

0.051 






0.081 

- 0.084 







1.743 


TABLE 10 

Optical depth between 50 GeV and 20 TeV eor sources at redshiet 
BETWEEN 0.01 AND 0 . 31 . 


TeV / z 

0.01 

0.03 

0.05 

0.07 

0.09 

0.11 

0.13 

0.15 

0.17 

0.19 

0.21 

0.23 

0.25 

0.27 

0.29 

0.31 

0.050 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.01 

0.01 

0.01 

0.01 

0.01 

0.061 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.01 

0.01 

0.01 

0.01 

0.01 

0.02 

0.02 

0.02 

0.03 

0.075 

0.00 

0.00 

0.00 

0.01 

0.01 

0.01 

0.01 

0.02 

0.02 

0.02 

0.03 

0.03 

0.04 

0.05 

0.05 

0.06 

0.091 

0.00 

0.00 

0.01 

0.01 

0.02 

0.02 

0.03 

0.04 

0.04 

0.05 

0.06 

0.07 

0.08 

0.09 

0.10 

0.12 

0.111 

0.00 

0.01 

0.02 

0.02 

0.03 

0.04 

0.05 

0.07 

0.08 

0.09 

0.11 

0.12 

0.14 

0.16 

0.18 

0.20 

0.136 

0.01 

0.02 

0.03 

0.04 

0.06 

0.07 

0.09 

0.11 

0.13 

0.15 

0.18 

0.20 

0.23 

0.26 

0.29 

0.32 

0.166 

0.01 

0.03 

0.05 

0.07 

0.09 

0.12 

0.14 

0.17 

0.21 

0.24 

0.27 

0.31 

0.35 

0.39 

0.44 

0.48 

0.202 

0.01 

0.04 

0.07 

0.11 

0.14 

0.18 

0.22 

0.26 

0.31 

0.36 

0.41 

0.46 

0.51 

0.57 

0.63 

0.69 

0.247 

0.02 

0.06 

0.11 

0.15 

0.21 

0.26 

0.32 

0.37 

0.44 

0.50 

0.57 

0.64 

0.71 

0.79 

0.87 

0.95 

0.302 

0.03 

0.09 

0.15 

0.21 

0.28 

0.35 

0.43 

0.51 

0.59 

0.68 

0.76 

0.86 

0.95 

1.05 

1.15 

1.25 

0.368 

0.04 

0.12 

0.20 

0.28 

0.37 

0.47 

0.57 

0.67 

0.77 

0.88 

0.99 

1.11 

1.23 

1.35 

1.48 

1.61 

0.450 

0.05 

0.15 

0.26 

0.37 

0.48 

0.60 

0.72 

0.85 

0.98 

1.11 

1.25 

1.40 

1.54 

1.69 

1.84 

2.00 

0.549 

0.06 

0.19 

0.32 

0.46 

0.60 

0.74 

0.90 

1.05 

1.21 

1.37 

1.54 

1.71 

1.88 

2.06 

2.24 

2.43 

0.671 

0.08 

0.23 

0.39 

0.56 

0.73 

0.90 

1.08 

1.27 

1.46 

1.65 

1.85 

2.05 

2.25 

2.46 

2.67 

2.88 

0.819 

0.09 

0.28 

0.47 

0.66 

0.86 

1.07 

1.28 

1.50 

1.72 

1.94 

2.17 

2.40 

2.63 

2.87 

3.11 

3.35 

1.00 

0.11 

0.32 

0.55 

0.77 

1.01 

1.24 

1.48 

1.73 

1.98 

2.23 

2.48 

2.73 

2.99 

3.25 

3.51 

3.77 

1.22 

0.12 

0.37 

0.62 

0.88 

1.14 

1.40 

1.66 

1.93 

2.20 

2.47 

2.74 

3.01 

3.28 

3.55 

3.82 

4.08 

1.49 

0.13 

0.41 

0.68 

0.96 

1.23 

1.51 

1.79 

2.07 

2.35 

2.63 

2.90 

3.18 

3.45 

3.72 

4.00 

4.27 

1.82 

0.14 

0.43 

0.71 

1.00 

1.28 

1.57 

1.85 

2.13 

2.41 

2.69 

2.97 

3.25 

3.52 

3.80 

4.07 

4.34 

2.22 

0.14 

0.43 

0.72 

1.01 

1.30 

1.58 

1.87 

2.15 

2.43 

2.72 

3.00 

3.29 

3.57 

3.86 

4.15 

4.44 

2.71 

0.15 

0.44 

0.73 

1.02 

1.31 

1.60 

1.90 

2.19 

2.49 

2.80 

3.11 

3.42 

3.73 

4.06 

4.38 

4.72 

3.31 

0.15 

0.45 

0.76 

1.07 

1.38 

1.70 

2.03 

2.36 

2.70 

3.05 

3.41 

3.77 

4.14 

4.52 

4.91 


4.05 

0.17 

0.50 

0.85 

1.20 

1.56 

1.93 

2.32 

2.71 

3.11 

3.52 

3.95 

4.38 

4.82 




4.94 

0.19 

0.59 

0.99 

1.41 

1.84 

2.29 

2.74 

3.21 

3.69 

4.18 

4.68 
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TABLE 10 — Continued 


TeV / z 

0.01 

0.03 

0.05 

0.07 

0.09 

0.11 

0.13 

0.15 

0.17 

0.19 0.21 0.23 0.25 0.27 0.29 0.31 

6.03 

0.23 

0.70 

1.18 

1.68 

2.19 

2.71 

3.25 

3.80 

4.36 

4.94 . 

7.37 

0.27 

0.82 

1.39 

1.98 

2.57 

3.19 

3.82 

4.46 



9.00 

0.32 

0.97 

1.63 

2.32 

3.03 

3.77 

4.53 




11.0 

0.38 

1.17 

2.00 

2.86 

3.77 

4.73 





13.4 

0.50 

1.56 

2.69 

3.89 







16.4 

0.73 

2.29 

3.98 








20.0 

1.11 

3.46 










TABLE 11 

Optical depth between 50 GeV and 6 TeV for sources at redshift 
BETWEEN 0.33 and 0.61. 


TeV / z 

0.33 

0.35 

0.37 

0.39 

0.41 

0.43 

0.45 

0.47 

0.49 

0.51 

0.53 

0.55 

0.57 

0.59 

0.61 

0.050 

0.01 

0.02 

0.02 

0.02 

0.02 

0.03 

0.03 

0.04 

0.04 

0.05 

0.05 

0.06 

0.06 

0.07 

0.08 

0.061 

0.03 

0.04 

0.04 

0.05 

0.06 

0.06 

0.07 

0.08 

0.09 

0.10 

0.11 

0.12 

0.13 

0.14 

0.15 

0.075 

0.07 

0.08 

0.09 

0.10 

0.11 

0.12 

0.13 

0.15 

0.16 

0.18 

0.19 

0.21 

0.23 

0.24 

0.26 

0.091 

0.13 

0.14 

0.16 

0.18 

0.20 

0.21 

0.23 

0.26 

0.28 

0.30 

0.33 

0.35 

0.38 

0.41 

0.43 

0.111 

0.22 

0.24 

0.27 

0.29 

0.32 

0.35 

0.38 

0.41 

0.45 

0.48 

0.52 

0.55 

0.59 

0.63 

0.68 

0.136 

0.35 

0.39 

0.42 

0.46 

0.50 

0.54 

0.59 

0.63 

0.68 

0.73 

0.78 

0.83 

0.88 

0.94 

0.99 

0.166 

0.53 

0.58 

0.63 

0.68 

0.74 

0.80 

0.85 

0.92 

0.98 

1.04 

1.11 

1.18 

1.25 

1.32 

1.40 

0.202 

0.76 

0.82 

0.89 

0.96 

1.04 

1.11 

1.19 

1.27 

1.35 

1.44 

1.52 

1.61 

1.70 

1.80 

1.89 

0.247 

1.03 

1.12 

1.21 

1.30 

1.40 

1.49 

1.59 

1.70 

1.80 

1.91 

2.02 

2.13 

2.24 

2.36 

2.47 

0.302 

1.36 

1.47 

1.58 

1.70 

1.82 

1.94 

2.06 

2.19 

2.32 

2.45 

2.58 

2.72 

2.86 

3.00 

3.14 

0.368 

1.74 

1.87 

2.01 

2.15 

2.30 

2.45 

2.60 

2.75 

2.90 

3.06 

3.22 

3.38 

3.54 

3.71 

3.88 

0.450 

2.16 

2.32 

2.49 

2.66 

2.83 

3.00 

3.18 

3.36 

3.54 

3.73 

3.91 

4.10 

4.29 

4.48 

4.68 

0.549 

2.62 

2.81 

3.00 

3.20 

3.40 

3.60 

3.81 

4.02 

4.23 

4.44 

4.65 

4.86 




0.671 

3.10 

3.32 

3.55 

3.77 

4.00 

4.22 

4.45 

4.69 

4.92 







0.819 

3.59 

3.83 

4.08 

4.32 

4.57 

4.82 










1.00 

4.03 

4.29 

4.54 

4.80 












1.22 

4.35 

4.62 

4.88 













1.49 

4.53 

4.80 














1.82 

4.62 

4.89 














2.22 

4.73 
















C. INDIVIDUAL SPECTRA 


The 106 VHE spectra studied in this paper are shown in Fig. [TBpO} The publications from which these data are 
extracted can be found in Table Observed VHE spectral points are shown in dark blue and the models best fitting 
these data are shown as a blue-gray solid line. The intrinsic spectral points are shown in light blue, after deabsorption 
with the best-fit EBL spectrum. The associated butterfly (one sigma envelope of the best-fit intrinsic model) is shown 
as a cyan contour. The butterflies associated with Fermi-hPCY data are shown in dashed gray for non-contemporaneous 
HE and VHE observations, and in dashed black whenever the HE and VHE data are quasi-contemporaneous. The 
border of each spectral panel is colored as in Fig. light blue indicating the gamma-ray cosmology sample, dark red 
the sources with under-constrained redshifts, and orange the FSRQs. 

The contemporaneity criterion is rather loose, and the HE and VHE integration windows differ in duration, albeit 
having a temporal overlap. For this reason, some VHE spectra show a flux level slightly larger than their quasi- 
contemporaneous counterpart (2 of 33) while some show a somewhat smaller flux level (4 of 33). No conclusion can 
be drawn from these mild discrepancies given the astrophysical uncertainties in the mechanism responsible for the 
variability of blazars and the differences between the HE and VHE exposures. 

T he constraints from the hardness of the intrinsic spectra based on quasi-contemporaneous HE observations (see 
Eq. 131 seem a posteriori rather weak. Nine spectra show non-zero contributions, the most important one being 
obtained for Mkn501_ARGD-YBJ_flare2011, with Ahe-vhe = 3.3 x 10“^. The remaining eight show values smaller 
than 10“^. 

The distribution of the intrinsic photon indices is shown in Fig. For the spectra described by log parabolas 
and exponential cut-off power laws, the indices are computed at the decorrelation energy. The minimum indices are 
F = 1.35 ±0.24 and F = 1.37 ±0.30, obtained for Mkn421_MAGIC_2006-04-27 and H1426+428_HEGRA_2002 respectively, 
indicating that all of these spectra are compatible with the maximum hardness, F > 1.5, expected within a standard 
synchrotron self-Compton model for an electron index of 2. 

All the spectra show good fits, with 103 of them having x^ probabilities larger than 15 %. The three poorest fits are 
PKS0301-243_HESS_2009-2010, Mkn501_ARG0-YBJ_flare2011, and Mkn421_VERITAS_2008JiighB with probabilities 
of 9%, 14%, and 14% respectively. 

Finally, we would like to stress a natural bias arising from the visual scanning of these spectra. It appears that 
some VHE spectra show important upward going fluctuations with respect to their best-fit model, see e.g. the last 
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Fig. 15. — Distribution of the intrinsic photon indices. 


lC310_MAGIC_2009-2010Jc 


IC310_UAGIC_2009-2010_higll-2=0.0189 


Mkn421_CAT_1998 • 2=0.031 


Mkn421_HEGRA_1999-2000 - 2=0.031 



Mkn421_HEQRA_2000-2001 -2=0.031 




Ukn421_T)BET_2000-2001 -2=0.031 




Ukn421_HESS_2004 - 2=0.031 





Mkn421_MAGIC_2006-04-22 - 2=0.031 


Mkn421_MAGIC_2006-04-24 - 2=0.031 


Mkn421_MAGIC_2006-04-25 - 2=0.031 


Mkn421_MAGIC_2006-04-26 - 2=0.031 



Mkn421_MAGIC_2006-04-27 - 2=0.031 




Mkn421_MAGIC_2006-04-28 - 2=0.031 







Fig. 16.— Best-fit spectra 1 of 4. See text for more details. 


spectral points of Mkn421_ARG0-YBJ_flux3 and lES0229-i-200_HESS_2005-2006. These cases are particularly striking 
to the eye with flux enhancements of ~ 1600 and ^ 23, respectively. In practice, such deviations are of rather small 
amplitude when normalized to the uncertainty on the flux, with deviations of only 1.2 cr and 1.3 cr for these two 
examples. The logarithmic scale of the plot is of course responsible for the visual bias that artificially amplifies upward 
going fluctuations. 
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Biteau & Williams 


Mkn421_ARGO-YBJ_flux2 • 2=0.031 


Mkn421_ARGO-VBJ_t!jx3 - 2=0.031 


Mkn421_ARGO-YBJ_flLx4 • 2=0.031 



Mkn421 _VERITAS_2008Jow - 2=0.031 



Ukn421_VERITAS_2008_tilghC • 2 = 0 .C 




Mkr501_ARGO-YBJ_2008-2011 - 2=0.034 




10 Iff'" 10"' 
Energy [TeV] 


1ES1959+650_HEGRA_2002 • 2=0.048 



BLLacertae_MAGIC_2005 - z=0.069 






Mkr501_HEGRA_1997 • 2=0.034 



Mkr501_ARGO-YBJ_nare20ll -2=0.034 




Energy [TeV] 

1ES1959-i-650_Whipple_2002 • 2=0.048 





Mkn421_VERITAS_2008_verylow - z=0.031 


Mkn421_VERITAS_200a_highA • 2=0.031 



Mkn421_TACTIC_2005-2006 • 2 = 



Mkn501_TACTIC_2005-2006 • 2=0.034 





Energy [TeV] 

1ES1959-^650_MAGIC_^006 - 2=0.048 



PKS2005-489_HESS_2004-2007 - 2=0.071 




Mkr421_VER(TAS_2008_highB - 2=0.031 



Mkn421_TACTIC_2008 - 2=0.031 



Mkn501_UAGIC_2006 - 2=0.034 



1 ES2344-t-514_MAGIC_2005-2006 - 2=0.044 




1 ES1959-i-650_VER;TAS_2007-2011 - 2=0.048 




Fig. 17.— Best-fit spectra 2 of 4. See text for more details. 
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SHBU00135S.9-18S40e_HESS_20M-2011 - z.0.095 



PKS2155-304_HESS_2002-06 • 2=0.116 



PKS2155-304_HESS_2003-08 • 2=0.116 



PKS2155-304_HESS_2006-07-0a -2=0.116 



PKS2155-304_HESS_2008-08-09 -2=0.116 




Energy [TeV] 


1RXSJ10101S.&-311909_HESS_2008-2010 -2.0.143 



RXJ0648.7-^1516_VERITAS_2010 - 2=0.179 



WComae_VER;TAS_2008-01-04 - 2=0.102 



PKS2155-304_HESS_2002-10 - 2=0.116 



PKS2155-304_HESS_2003-08 - 2=0.116 



PKS2155-304_MAGIC_2006-07_08 - 2=0.116 



B32247t381_MAGIC_2010 - 2=0.119 



1 ES0806-t-524_VERITAS_2006-2008 - 2=0.138 



H2356-309_HESS_2004 -2=0.165 



1 ESI 218 t304_VERITAS_2007 - 2=0.182 



1ES1312-423_HESS_2004-2010 - 2=0.105 



PKS2155-304_HESS_2003-06 - 2=0.116 



PKS2155-304_HESS_2003-10-11 - 2=0.116 



PKS2155-304_HESS_2006-08-02 -2=0.116 



RGBJ0710-^591_VERITAS_200a-2009-2=0.125 



1 ES0229-f200_HESS_2005-2006 - 2=0.14 





VERJ0521+21 1_VER;TAS_2009-2010 - 2=0.108 



PKS2155-304_HESS_2003-07 - 2=0.116 



PKS2155-304_HESS_2005-2007 - 2=0.116 



PKS2155-304_HESS_2006-0a-03 - 2=0.116 



H1426t428_HEGRA_1999_2Q00 - 2=0.129 



Energy [TeV] 
H2356-309_HESS_2006 -2=0.165 



1 ESI 101-232_HESS_2004-2005 - 2=0.186 



Fig. 18. 


Best-fit spectra 3 of 4. See text for more details. 




























































































Energy flux vF„ [erg ernes’] Energy flux vF„ [ergem'^s'] Energy flux vF„ [ergem'^s'] Energy flux vF„ [ergem'^s’] Energy flux vF„ [erg ernes'] Energy flux vF„ [ergem'^s^ Energy flux vF„ [erg ernes’] 


28 
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1ES0347-121_HESS_2006-09-12 • 2=0.188 



1ES1215+303_MAGIC_2011 -2=0.237 







1 ES1215-i-303_VER;TAS_2008-201 2 • 2=0.237 




1 ES0414+009_VERITAS_2008-2011 - 2=0.287 




3C66A_MAGIC_2009-2010 • 2=0.335 



PKS0447-439_HESS_2009-2010 • 2=0.5 



PKS1510-089_HESS_2009 - 2=0.361 



PG1553t113_HESS_2005-2006 • 2=0.433 



PG1553-fl 13_MAGIC_2007 • 2=0.^ 



PG1553-t-113_MAGIC_2008 - 2=0.4 



PGl553-f113_MAQIC_2009 • 2=0.433 



PKSl424t240_VERITAS_2009 - 2=0.604 



PG1553-^113_VERITAS_2010•2012-2=0.433 



PKS1424-^240_VERITAS_2011 -2=0.604 



PG1553-t-113_HESS_2012_nare - 2=0.433 




3C279_MAGIC_2006 - 2=0.5362 



PKS1424+240_MAGIC_2009 - 2=0.604 



PKS1424+240_MAGIC_2010-2=0.604 PKS1424+240_MAGIC_2011 -2=0.604 




Fig. 19.— Best-fit spectra 4 of 4. See text for more details. 























































